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HI-ALPHA FOREBODY DESIGN: PART I 
METHODOLOGY BASE AND INITIAL PARAMETRICS 

William H. Mason 
R. Ravi 

Virginia Polytechnic Institute and State University 
SUMMARY 

The use of Computational Fluid Dynamics (CFD) has been investigated for the analysis and 
design of aircraft forebodies at high angle of attack combined with sideslip. The results of the 
investigation show that CFD has reached a level of development where computational methods can 
be used for high angle of attack aerodynamic design. The classic wind tunnel experiment for the 
F-5A forebody directional stability has been reproduced computationally over an angle of attack 
range from 10° to 45°, and good agreement with experimental data was obtained. Computations 
have also been made at combined angle of attack and sideslip over a chine forebody, demonstrating 
the qualitative features of the flow, although not producing good agreement with measured 
experimental pressure distributions. The computations were performed using the code known as 
cf!3d for both the Euler equations and the Reynolds equations using a form of the Baldwin- 
Lomax turbulence model. To study the relation between forebody shape and directional stability 
characteristics a generic parametric forebody model has been defined which provides a simple 
analytic math model with flexibility to capture the key shape characteristics of the entire range of 
forebodies of interest, including chines. 
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INTRODUCTION 


High angle of attack aerodynamic characteristics are key to agility of advanced fighter 
aircraft. Considerable attention has been devoted to flow asymmetries, advanced control device 
concepts (forebody ’’flaps”), and modem control system design (ref. 1-4). However, the 
requirement for advanced control devices and control systems is in large part due to directional 
stability characteristics. If the basic directional stability characteristics are aerodynamically tailored 
by proper shaping of the design, demands on control systems and control devices can be 
dramatically decreased, and complicated forebody mechanisms (which are undesirable from an 
aircraft designer's point of view) may not be required. 

The aerodynamic design problem can be described using typical aerodynamic characteristics 
of fighter aircraft. Figure 1 illustrates the problem in terms of typical characteristics of an 
advanced fighter. At low angle of attack the vertical tail provides directional stability. As the angle 
of attack increases the tail loses effectiveness because it is located in the separated flow region 
behind the fuselage and wing. For some configurations the directional stability increases again as 
the angle of attack increases, as shown in the figure. Figure 2 is taken from the report by Grafton, 
et. al., (ref. 5). In this case the directional stability of the forebody alone is nearly identical to the 
directional stability of the entire aircraft above 30° a, and shows the dominant role played by the 
forebody in causing the directional stability to start increasing for an F-5A type of configuration. 
Figure 3 (ref. 6) provides a schematic of a typical forebody and flowfield responsible for favorable 
characteristics, although in this case the concept is quantitatively in error. The leeward vortex is 
actually bigger but farther away from the surface than the windward vortex. The low pressure 
associated with the vortex located very close to the surface acts to “pull” the body back to a smaller 
sideslip, and thus provides a stabilizing moment. The possibility of controlling the flow is shown 
in the figure, where the cross section is not axisymmetric, and has a small radius of curvature or 
crease at the max half breadth line location. Figure 4 shows the desirable changes to the C n ^ 
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characteristics shown in figure 1. Ideally the loss of stability should be delayed and the rapid 
variation should be minimized, with a resulting maximum C n ^ being less than the extremely high 
values sometimes observed. High values of C n ^ are frequently considered to be an indication of 
poor C Kr characteristics. 

Previous development work in this area has been carried out experimentally. Figure 5 
illustrates a range of forebody shapes that have been considered for use on fighter aircraft. 
However the capability now exists to use computations to investigate the design of aircraft 
components at high angle of attack. This has been demonstrated by recent computational results. 
One practical example is the analysis of the F/A-18 (ref. 1 1). Detailed examinations have also been 
made for tangent ogives, including the effects of surface perturbations on flow asymmetry (ref. 
12-14). However, none of the previous work has examined the capability of the methods to 
compute sideforces arising from combined angle of attack and sideslip on non-axisymmetric 
forebodies. 


Although the general problem of obtaining positive directional stability levels for fighter 
aircraft at high angle of attack is complicated and component interactions are important, the 
forebody is crucial, as shown in figure 2, and because of its location necessarily less affected by 
interactions than other components . Therefore, the broad objective of the current effort is directed 
toward developing an understanding of forebody static directional stability aerodynamics and 
methods to tailor the isolated forebody characteristics for desirable high angle of attack 
characteristics. Specifically: i) how can a favorable contribution to directional stability be induced at 
lower angles of attack and ii) how can extremely large values of C n p (which presumably are 
indicative of poor C n ^_ characteristics) be eliminated? 

This report addresses the necessary first steps required to achieve the overall objective: the 
capability and credibility of using an advanced computational code to conduct initial aerodynamic 
design of forebodies. To investigate the current capabilities of CFD, a current, widely used code, 
cf!3d (ref. 15), was used to compute two cases for which experimental data is available: the F-5A 
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forebody tested by Grafton, et. al. (ref. 5), and the chine forebody tested by Erickson and Brandon 
(ref. 7). A discussion of the results is preceded by a description of the initial work done to use 
cf!3d in a combined a / ft flowfield. After the comparison with the experimental results, a 
parametric model for a generic forebody is presented which is completely defined by simple 
analytic formulas. This generic forebody model is capable of describing forebody geometries over 
an extremely large class of shapes, from sharp edge chines to rectangular cross sections. 
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LIST OF SYMBOLS 


a 

b 


b' 



m,n 



Rei 

s ref 

* 

u 


x,y,z 


maximum half breadth of the generic forebody definition 
maximum centerline distance of the generic forebody definition 
wingspan 

mean aerodynamic chord 
lift-force coefficient, lift 

pitching-moment coefficient, pitching moment t/q 0 J> re jc 
yawing-moment coefficient, yawing moment Jq^^jb' 
directional stability derivative, dC n / d(5 
pressure coefficient, (p-Poo)fQoo 
axial-force coefficient, axial fore e/q^^ 
side-force coefficient, side force/q^^ 
normal-force coefficient, normal force/g^^y 
local side-force, section side foTce/q c JS re j- 
diameter of circular cross section of tangent ogive 
fuselage station 

adjustable parametric coefficients 
free stream Mach number 
model length 

Reynolds number based on model length, / 
reference area 

wall friction velocity, Vt w / p 

body coordinate system : x positive aft along model axis, 
y positive to right and z positive up 
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moment reference center 

jjc 

inner law variable, yu /v 
angle of attack, deg 
angle of sideslip, deg 

azimuthal angle, measured clockwise from windward plane at any cross 
section 


COMPUTATIONAL BASIS 


The baseline code for this work was the NASA Langley program cf!3d. At the start of the 
effort version 1.1 of this code was selected as being appropriate for this investigation. In this code 
the three dimensional compressible viscous flow around the body is found by solving the 
conservative form of the dimensionless, thin layer Navier-Stokes equations. 

For turbulent flow cfI3d uses the Reynolds averaged counterparts of the Navier-Stokes 
equations. An algebraic turbulence model is used, where the turbulent viscosity fi t is obtained by 
using the two layer algebraic eddy viscosity of Baldwin and Lomax (ref. 16) as modified by 
Degani and Schiff to account for the special characteristics of strongly separated flows in 
references 17 - 19. The Roe type inviscid flux difference splitting scheme option in the code was 
used to obtain the results presented in this report 

The farfield computational boundaries were chosen sufficiently far away from the body so as 
not to affect the forces on the body. They were chosen to be consistent with those used by other 
investigators e.g.,F/A-18 (ref. 11) and tangent ogive (ref. 12). The exact distances have been 
given under each specific case discussed later. Inflow-outflow boundary conditions were used at 
the farfield inflow, farfield outflow and the farfield outer boundary. For turbulent computations the 
grid ahead of the nose was treated as a separate block (as described later). When viscous effects are 
included the no-slip as well as non-penetration conditions are enforced on the body. At the plane of 
geometric symmetry periodic conditions are used to include sideslip as well as angle of attack in the 
freestream flow. The temperature boundary condition is treated by defining the body to be 
adiabatic. 

The computations were made on the NASA Langley Cray-2 computer. Some initial startup 
activity was required before beginning the forebody study. The available version of the code was 
modified, as described in Appendix A, to handle the boundary conditions for combined a/p 
flows, and then validated for operation at combined angles of alpha and beta. 
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The primary, and minor, modification of the code was the incorporation of a new boundary 
condition. For the work conducted here a full grid was used wrapping all the way around the 
body, instead of a typical plane of symmetry grid. Therefore, instead of plane of symmetry 
boundary conditions, the grid interface in the cross flow plane was setup to allow flow across the 
grid boundary on the leeward side. 

The computational grids were generated using stacked 2D grids at each axial location. The 2D 
grid in the cross plane was generated using a grid generator provided by W. McCrory (ref. 20) that 
was originally developed to construct the grid for the calculation of the SR-71 flowfield (ref. 21). 

Analysis of flowfield calculations was carried out for the more complicated aspects of the 
grid and flow visualization using the NASA Ames code PLOT3D (ref. 22). This program was 
used on the Aerospace and Ocean Engineering Department Iris 4D graphics workstation. 


LAMONT TANGENT OGIVE 

To check the code and grid setup for this work a simple geometry extensively studied by 
other investigators was selected. The so-called “Lamont tangent ogive” (ref. 23) was used for the 
work. This is a 3.5 diameter tangent ogive nose with a cylindrical afterbody. The total length of the 
ogive and cylindrical forebody was 25 diameters. The objective of computing the flow over this 
relatively simple geometry was to test the minor modifications that were made to cfl3d, including 
the added boundary condition. The ability of the code to handle combined angle of attack and 
sideslip is a key part of the study and needed to be verified with a geometry which would provide 
the same results for a at zero p and p at zero a. Inviscid computations were made at 

M„ = 0.2 for two different cases, one for a = 20° and no sideslip and the other with p = 20° with 

» 

no angle of attack. 
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Tangent Ogive Forebody Grid 


The three dimensional grid was constructed from two dimensional O-type cross flow grids 
which are longitudinally stacked, constituting an H-0 topology. The boundaries upstream of the 
nose and radially away from the body were selected to be consistent with similar calculations done 
by Hartwich and Hall (ref. 12). The forward boundary extends upstream of the nose by 0.502/ 
and radial outer boundary extends 1.06/ from the model centerline. The reference length / is 
equal to the longitudinal extent of the forebody, which was 25 diameters. 

The grid used for the initial inviscid calculations had 24 points in the radial direction and 57 
points in the full circumferential direction. Longitudinally, the grid was clustered near the nose 
with 22 stations on the forebody and 9 stations ahead of the nose. The grid upstream of the nose 
was also longitudinally stretched to provide resolution near the nose. The entire inviscid grid 
consisted of 42,408 points. The results obtained using this grid resulted in irregularities in the 
surface pressures at the ogive-cylinder intersection and downstream to be discussed below. The 
problem was reduced by increasing the grid resolution. Thus, subsequent calculations were made 
by increasing the number of grid points in the circumferential direction to 69 points with total 
number of points in the entire grid increasing to 51,336 points. 

Results and Discussion of the Lamont Tangent Ogive Computations 

Initial inviscid calculations used the 24 x 57 x 31 grid for the case of a = 20° and P = 0°. 
These results are shown in figures 6 and 7. The experimental data shown in these figures are due 
to Lamont (ref. 23). Here they are taken from the paper by Hall and Hartwich (ref. 12) and include 
an adjustment in the pressures. The data were obtained at a wind tunnel freestream Reynolds 
number of 0.2 x 10 6 based on the diameter of the tangent ogive. Figure 6 shows the surface 
pressures on the windward plane of symmetry. The difference in the surface pressure between the 
experimental and the computational data at sections close to the nose is attributed to the poor 
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resolution of the stacked O-grid type topology at the nose. Figures 7(a) to 7(d) show the surface 
pressure coefficient at the axial stations of x/d = 0.5, 2.0, 3.5 and 6.0. The first two locations are 
on the ogive forebody, the third is at the junction between the forebody and the cylindrical 
afterbody and the fourth is on the cylindrical afterbody. The 6 shown on the x axis is the angle 
measured clockwise from the windward plane at any cross section with 0=0° corresponding to 
the windward plane and 6 = 180° corresponding to the leeward plane. The inviscid computations 
agree fairly well with the data on the windward side. Agreement is poor near the nose (fig. 7(a)). 
Presumably this is due to a lack of grid resolution, as seen in figure 6. The agreement is 
surprisingly good at x/d = 2.0 considering the expectation that viscous effects would have 
significantly altered the pressure distribution at this angle of attack. However, agreement does 
deteriorate further downstream, where viscous effects do begin to make a significant effect on the 
surface pressure distribution as shown in figures 7(c) and 7(d). The irregularities in the computed 
surface pressures near 6 = 90° decreased significantly when the 24 x 69 x 31 grid was used. All 
the following figures incorporate this finer grid. 

Using the refined mesh, computations were also made for the case of a = 0° and /? = 20° to 
verify that the surface pressures and forces were equivalent to the a = 20° and ft = 0° case. The 
surface pressures at different cross sections are shown in figures 8(a) to 8(c) for both the cases. 
The ordinates for the case of a = 0° and P = 20° have been shifted by 90° for comparison to the 
surface pressures from the a = 20° and P = 0° case. The pressures are identical in both cases. 
Table [ 1 ] gives the forces for these cases. The forces and moments are nearly identical with the 
orientation of the body in each case. Only small differences occur between components. The 
replacement of the plane of symmetry boundary condition with the periodic boundary condition 
also weakened the convergence. 

This computation establishes the validity of the new boundary condition and the ability to 
treat side slip with angle of attack. Having established confidence in the code for sideslip as well as 
angle of attack calculations, the calculation of aerodynamic interest were initiated. 
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F-5A FOREBODY 


The flowfield around the F-5A forebody was computed to study the use of an advanced code 
like cfI3d in predicting the directional stability characteristics of forebodies at high angles of 
attack. This forebody had been tested by Sue Grafton , et.al. at NASA Langley Research Center 
and the results are available in ref 5. 

F-5A Forebody Geometry Math Model 

The computational math model was constructed using classical conic lofting techniques (ref. 
24) based on the lofting information provided by Tom Heglund of Northrop. Figures 9(a) to 9(e) 
show the F-5A forebody geometry with cross section shapes at different locations along the length 
of the body, together with the data required to construct the numerical model. The configuration 
was modeled from the nose longitudinally back to x = 225 in. The computational model did not 
include the canopy. At sections close to the nose, as shown by figures 9(b) and 9(c), the upper and 
the lower maximum half breadth points coincide. However, at sections downstream the upper and 
lower maximum half breadth points are separated by a flat side wall. Because of this, grid points 
were clustered near the maximum half breadth points forward of the flat side wall to provide 
adequate definition on the flat wall portion of the forebody and create smoothly varying 
longitudinal grid lines. The computer code to generate the F-5A surface co-ordinates is given in 
Appendix B. 

To verify that the wind tunnel model geometry agreed with the math model, Sue Grafton of 
NASA Langley provided the forebody from the wind tunnel model of the F-5A. The model was a 
0.17-scale model of the full scale airplane. Figures 10(a) and 10(b) show the wind tunnel model 
forebody and the setup that was used to make templates of the cross sections. The measured shape 
of the forebody at different cross sections was then digitized from a tracing of the shape defined by 
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the template and compared with the computational model. The cross sectional shapes of the 
computational and the wind tunnel model agree very well considering the method used to obtain the 
model shape, as shown in figures 11(a) to 11(f). The shaded surface pattern of the F-5A 
computational model is given in figure 12. 

F-5A Grid 

The three dimensional grid was constructed from two dimensional O-type cross flow grids 
which are longitudinally stacked, constituting an H-0 topology. The boundaries upstream of the 
nose and radially away from the body were initially selected to be consistent with the F/A-18 
calculations in ref 11. These were later found to be sufficiently far away from the body so as not to 
affect the forces and moments on the body. The forward boundary extends upstream of the nose 
by 0.62/ and radial outer boundary extends 0.98/ from the model centerline. The downstream 
boundary, as explained earlier, extended from the nose back to x = 225 in ( in full scale ) which 
is as shown in figure 9(a) without the canopy. The reference length / is equal to the longitudinal 
extent of the forebody which was 31.025 in using the wind tunnel model scale. 

For inviscid calculations the grid had 45 points in the radial direction and 93 points in the full 
circumferential direction. This number was increased based on the experience with the tangent 
ogive and the need to resolve the details of the cross section. Longitudinally, the grid was clustered 
near the nose with 25 stations on the forebody and 8 stations ahead of the nose. The grid upstream 
of the nose was also longitudinally stretched to provide resolution near the nose. The entire 
inviscid grid consisted of 138,105 points. 

For viscous calculations the grid had 65 points in the radial direction with longitudinal and 
circumferential grid points remaining the same as used for the inviscid calculations. The total 
number of grid points for the viscous grid was 199,485 points. The grid upstream of the nose was 
treated as a separate block for turbulent flow calculations. Figures 13 to 15 show the grid used for 
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the viscous calculations at different cross sections downstream from the nose. The reason for the 
clustering of the grid points near the maximum half breadth point at stations close to the nose was 
described above. Figure 16 shows the longitudinal stretching of the grids. For the grid that was 
used, it can be seen that good resolution of the grid at the nose is required to resolve the flow 
details there.The actual grid used in the computations had more stations near the nose than that 
shown in figure 16. 

The grid was established with sufficient normal clustering near the surface to adequately 
resolve the laminar sublayer for the turbulent boundary layer flow. This grid produced an average 
normal cell size of approximately 10 " 4 /. At the wind tunnel ffeestream conditions ( = 0.2, 
Rei = 1.25 x 10 6 , and a = 20 ) the baseline grid typically resulted in a value of y + ~ 2 at the 
first mesh point away from the body. 

Results and Discussion of Computations on the F-5A Forebody 

Inviscid calculations were performed for a = 20°, 30°, 40° and 50° with side slip angles of 
p = 0°, 5° and 10°. Turbulent calculations were performed for a = 10°, 20°, 30°, 40° and 45° 
and P =5°. All computations were performed on the Cray-2. The inviscid computations took 
approximately 0.55 hours of CPU time for convergence, and the turbulent computations took 
approximately 4.5 hours to reduce oscillations in the lift coefficient C £ and the yawing moment 
coefficient C n to a negligible level. 

Convergence histories for a typical inviscid case is given in figure 17 in terms of residual and 
lift. The convergence histories for a typical viscous computation are shown in figure 18 which 
includes the residual, the lift force, side force and the yawing moment coefficients. The inviscid 
calculations which were started with a CFL number of 2 and ramped over 100 steps to a CFL 
number of 6 demonstrated strong convergence. Unlike the inviscid case, the turbulent 
computations show low frequency decaying oscillations in values of Cj^ and C n and require 
about 5000 iterations for a fully converged solution. The turbulent computations were started with 
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a CFL number of 0. 1 and ramped over 100 steps to a CFL number of 0.3 The detailed behavior of 
the residual and forces change abruptly after 3000 iterations as shown in figures 18(a) and 18(b). 
This is the result of reducing the CFL number once the solution ceased to converge, thus, making 
the oscillations decay further to a steady solution. 

The F-5A forebody directional stability experimental data from ref 5 are shown along with 
the computed inviscid and viscous results in figure 19. The computed values of C np were 
obtained as the finite difference of the yawing moment coefficients C n . These values being small, 
it was required that the solution be fully converged to avoid errors in the computed values of C n p 
It can be seen that both inviscid as well as viscous results show the stabilizing effect of the 
forebody at high angles of attack. At low angle of attack (10° and 20°) the inviscid and viscous 
calculations agree with each other and the data. At 30° the viscous computations continue to agree 
with the data, and at the higher angles of attack (40° and 45°) differences begin to appear. This is 
reasonable considering both the thin layer assumption and the accuracy of the turbulence model for 
massively separated flows. Nonetheless, the ability of the computations to reproduce the 
experimentally observed contribution of the F-5A forebody to directional stability is an exciting 
result. The viscous calculation results establish the effectiveness of an advanced code like cfI3d in 
computing the directional stability characteristics for forebodies at high angles of attack. However 
the inviscid results are in error. The Euler solution results are shown below to contain the spurious 
non-physical flow separation arising due to lack of grid resolution. The grid resolution effects were 
not studied in this report but the issue should be considered as a part of future study. 

Figures 20 and 21 show the inviscid and turbulent surface pressure distributions at various 
cross section stations, respectively. The corresponding cross sectional shapes are shown at the 
bottom of the figure. The results presented here were carried out on the respective standard grids 
for the case of a = 40° and (5 = 5°. The negative peak pressures are due to the vortices on the 
upper surface of the cross section and are shown more clearly in the surface flow visualization 
pictures presented in the following paragraph. The asymmetry in the pressure distribution due to 
the side-slip can also be seen. Figure 22 provides a comparison of the inviscid and turbulent 
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results. At station 14 the viscous results clearly show the effect of the vortices, with two low 
pressure regions underneath the windward and leeward primary vortices. Curiously, at the aft 
station (29.61), the inviscid results contain the low pressure peaks. At this station the true vortices 
in the turbulent calculations are well away from the surface, and the flow is massively separated, 
with the associated low pressure gradients along the surface. The three low pressure peaks in the 
turbulent solution are due to the windward primary vortex, the leeward primary vortex and the 
leeward secondary vortex. The inviscid results contain low pressure regions that are related to the 
cross section curvature. 

The flowfield results shown in figures 23 to 27 have been generated using computational 
flow visualization techniques (PLOT3D). Figure 22a shows the inviscid surface flow pattern for 
a = 40° and {5 = 5°. The dye is injected at two axial stations to observe the formation of 
separation lines, if any. The inviscid separation lines seem to indicate some kind of spurious 
separation occurring on the surface. The turbulent flow surface oil flow pattern, is shown in figure 
23b for the same a and j3. Here, primary separation lines starting from near the nose as well as 
the secondary separation lines starting from further downstream are shown. 

Figures 24 and 25 show the cross sectional pressure contours at the axial station x = 14.025 
in. from the nose for the inviscid and viscous flows respectively, where these results are again for 
a = 40° and /J = 5°. The magnitudes associated with the contour quantities are displayed with a 
color bar on the left.The viscous leeward vortex is bigger and farther away from the surface than 
the windward vortex in figure 25. The lower surface pressure associated with the vortex closer to 
the surface acts to pull the body back to smaller sideslip, and thus provides a stabilizing moment. 
The vortex pattern and asymmetry are more clearly seen in the figures 26 and 27 which are shaded 
based on the stagnation pressure. The leeward vortex in the case of turbulent flow is farther away 
from the surface than that of the inviscid flow. This difference in the relative positions of the 
vortices in turbulent and inviscid cases is the reason for the different directional stability shown in 
figure 19. The turbulent flow also shows the secondary vortices being shed from under the 
primary vortices further downstream from the nose. 
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ERICKSON CHINE FOREBODY 


The surface pressure distributions from the chine forebody tested by Erickson and Brandon 
(ref. 7) were used to investigate the capability of CFD to compute the flow over a chined forebody. 
The experimental configuration is shown in figure 28 (ref. 7). The geometry information used to 
construct the grid was supplied by Robert Hall of NASA Langley, and the pressure data was 
supplied by Gary Erickson, also of NASA Langley. 

The wing tunnel model had rows of pressure orifices at FS7.19, FS 13.56, and FS 19.94. All 
of the data was taken with the wing present. It was assumed that the presence of the wing would 
not affect the results at the first station, and have minor effects at the second station. The third row 
of pressures was considered too close to the wing to provide valid comparisons with the forebody 
calculations, but are presented to provide additional information. 

Erickson Chine Forebody Math Model 

The Erickson forebody definition was not available in a convenient form for development of 
the surface definition. A CAD drawing was made available, together with digitized ordinates of the 
cross-sections contained on the drawing. We were also supplied with digitized coordinates of the 
drawing generated by Nielsen Engineering and Research, Inc. As a basis for the shape, Robert 
Hall supplied a copy of his notes defining the class of theoretical shapes intended to be used to 
design the forebody. The forebody was required to blend smoothly to an existing NASA Langley 
model, and this requirement led to deviations from the theoretical concept 

The surface is constructed using standard lofting techniques. First, the key longitudinal lines 
are defined as a function of the axial station. In this case the lines are the top center line, the bottom 

* The CAD drawing of the model forebody was identified as “Roll Stability Model 
Forebody”, 11 July 1984, C.J. Rozo, Change A 
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center line, and the max half breadth line. A cross section model is then developed which connects 
these lines smoothly and produces a section which varies smoothly with axial location. For this 
model the upper and lower cross sections are developed independently of each other. The listing of 
the FORTRAN code of the math model is given in Appendix C. 

Longitudinal Lines 

Upper surface centerline: from the nose to FS7.5 the top centerline is defined by a portion of 
a 19.5° tangent ogive. Aft of this station the shape is defined by combination of three cubic spline 
segments and two straight line segments. This model was based on careful analysis of the CAD 
drawing. 

Lower surface centerline: from the nose to FS6.0 the bottom centerline was defined by a 
portion of a 15° tangent ogive. Aft of this station the shape was defined by a combination of two 
cubic line segments and two straight line sections. Again, the model was developed based on an 
analysis of the CAD drawing. 

Max half breadth line: from the nose to FS8.0 the max half breadth line was defined by a 
portion of a 27.5° tangent ogive. A cubic spline was used from FS8.0 to FS19.0, followed by a 
zero slope straight line segment. 

In all cases the line segments and spline end point conditions were chosen to provide 
smoothly varying changes in the slope of the lines, and to minimize jumps in the curvature 
distribution. Analysis of the CAD drawing suggests that the modifications to the theoretical shape 
originally outlined by Hall resulted in some abrupt changes in curvature distribution in the actual 
model lines (recall that a point and slope type match between surfaces will normally result in a 
discontinuity in curvature at that point). The exact details of the models are easily found in the 
listing in Appendix C. 
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Goss Sections 


Upper surface cross section: Forward of FS 7.19 the Hall equation was found to fit the 
supplied digitized ordinates from the drawing well. Aft of the FS7.17 the section was modified to 
match the existing circular fuselage at FS23. Thus another station was defined at FS23, which was 
initially circular, and then departed to make the chine using cubic splines. Between these two 
sections a combination of the two function was used based on the blending function method 
described by Barger and Adams (ref.25) 

Lower surface cross section: When normalized by the max half breadth line and the lower 
centerline the digitized representation of the cross sections was found to be the same at all cross 
sections. Thus only one cross sectional shape was required. This cross section was constructed 
using Hall’s equation as the basis, and adding a spline curve defining the difference between Hall’s 
equation and the actual shape. The difference was found by subtracting the analytic shape given by 
Hall from the measured ordinates. Because the measured coordinates were slightly noisy, a smooth 
curve was developed by selecting a small number of points from the data, and using a spline 
interpolation between the points to evaluate intermediate values. 

The results of the math model are compared with the ordinates measured from the drawing in 
figure 29 for the full range of stations. Figure 30 presents details at the three stations where 
pressure data was obtained. This is the best smooth model of the forebody we were able to 
construct. The surface shaded PLOT3D view of the model is shown in figure 31, and demonstrates 
the resulting smooth contour of the math model. Construction of an accurate math model for the 
Erickson forebody was one of the most time consuming parts of the work described in this report. 
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Erickson Chine Forebody Grid 


Using a procedure similar to the approach used for the F-5A grid, the three dimensional grid 
around the Erickson forebody was constructed from two dimensional O-type cross flow grids 
which were longitudinally stacked, constituting an H-0 topology. The boundaries upstream of the 
nose and radially away from the body were again selected to be consistent with the F/A-18 
calculations in ref. 11. The downstream boundary was extended from 19.94 in to 30.00 in with 
the same cross section. The forward boundary extends upstream of the nose by 0.61/ and radial 
outer boundary extends .99/ from the model centerline. The reference length / is equal to the 
longitudinal extent of the forebody (30.00 in). 

The inviscid calculation grid had 45 points in the radial direction and 101 points in the full 
circumferential direction. Longitudinally, the grid was clustered near the nose with 25 stations on 
the forebody and 8 stations ahead of the nose. The grid upstream of the nose was also 
longitudinally stretched to provide resolution near the nose. Care was taken to have axial stations at 
sections where the experimental surface pressures were available. These were at a distance of 7.19, 
13.56 and 19.94 inches from the nose along the length of the body. The entire inviscid grid 
consisted of 149,985 points. 

As compared to the inviscid grid, the viscous calculation grid had 65 points in the radial 
direction with longitudinal and circumferential grid points remaining the same. The total number of 
grid points for the viscous grid was 216,645 points. As described for the F-5A calculations, the 
grid upstream of the nose was treated as a separate block for turbulent flow calculations. Figure 32 
shows the viscous grid structure at a typical cross section. 

The baseline grid was established with sufficient normal clustering near the surface to 
adequately resolve the laminar sublayer for the turbulent boundary layer flow. This grid produced 
an average normal cell size of approximately 10'^/. At the wind tunnel freestream 
conditions ( M ^ = 0.2, Rei = 1.02 x 10 6 , and a = 20° ) the baseline grid typically resulted in 
a value of y + » 2 at the first mesh point away from the surface. 
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Results and Discussion of the Computations on the Erickson Chine Forebody 


Inviscid and turbulent calculations were initially performed for a =30° with side slip angles 
P = 0° and 10°. Inviscid computations took approximately 0.75 hours of CPU time for 
convergence to five orders of magnitude.The turbulent computations took approximately 4.0 CPU 
hours to reduce oscillations in lift force coefficient C L and yawing moment coefficient C n to a 
negligible level, by which time the residual went down two orders in magnitude. 

The convergence summary for a typical inviscid and viscous case are given in figures 33 and 
34 respectively. Convergence histories of the residual and the normal force coefficients are shown. 
As in the case of F-5A the inviscid calculations were started with a CFL number of 2 and ramped 
over 100 steps to a CFL number of 6. Unlike the inviscid case, the turbulent computations show 
low frequency decaying oscillations in values of C L and require about 4000 iterations for a fully 
converged solution. The turbulent computations were started with a CFL number of 0.1 and 
ramped over 100 steps to a CFL number of 0.3. Once the solution ceased to converge after 2000 
iterations the CFL number was increased and the oscillations decayed further to a steady solution. 
This change in CFL number is the reason for the sudden change in the behavior of the plot in 
figure 34. 

Figures 35 and 36 present the computed inviscid surface pressure distributions on the upper 
surface at the three pressure instrumented cross sections, with the cross sectional shapes shown 
below each of the plots . The experimental data were obtained from Gary Erickson. The details of 
the experimental investigation are available in ref 7. The axial stations where the pressures have 
been plotted were at a distance of 7.19, 13.56 and 19.94 inches from the nose along the length of 
the body. Figure 35 shows the pressures for the a = 30° and (5 = 0° case. At the section closest 
to the nose (7.19 in.) the computations predict the pressure level very well near the center section 
of the body, and underpredict the surface pressure suction levels outboard, where the suction peak 
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occurs. Moving downstream, the trends are the same, except that the agreement at the center 
section also becomes less accurate. Figure 36 shows the surface pressure comparisons for the case 
of a = 30° and ft = 10°. In this case the windward pressure suction levels, as well as the center 
section pressures are predicted very well by the inviscid computations. The agreement deteriorates 
at the FS 19.94 station, where wing effects not included in the calculation would become important. 
The leeward pressure levels are poorly predicted at all the stations. 

The effect of viscosity on the surface pressure comparisons presented above are given in 
figures 37 and 38. Including viscosity effects does not improve the agreement. Figure 37 provides 
the results at zero sideslip. Viscous effects did not change the pressure levels at the mid section of 
the forebody, but did have a large effect on the peak suction pressure level, reducing the peak 
pressures and resulting in even worse agreement with the experimental data. The same trend is 
observed for the results at 10° sideslip, as shown in figure 38. 

Due to the disappointing agreement between the experiment and computations, an estimation 
of the possible effects of the wing on the forebody flowfield was made. To investigate the effect of 
the wing on the flow over the forebody Lamar's Vortex Lattice code ( VLM4997, ref 26 ) was run 
for an isolated wing model of the wind tunnel model for the case of a = 30° and p = 0°.The 
velocity field due to the wing lift was then computed at the pressure measurement locations. 
Although approximate, this calculation provides an estimate of the additional induced angle of 
attack that would exist on the forebody due to the presence of the wing. Two models of the isolated 
wing were considered, as shown in figure 39. The first model had the leading edge sweep 
extended to the centerline and the second model had the wing nose clipped at the wing-chine 
forebody intersection. Using the experimental value of the lift at a = 30°, C L = 1.0, the code 
computed an angle of attack of 26°. Under these conditions the induced angles of attack as 
computed by the VLM code are given in figure 40. The effect is approximately 1° at FS7.19, 2° at 
FS 13.56, and becomes large at FS 19.94, which is very close to the wing. 

Based on the induced angle of attack analysis, an inviscid solution was computed using an 
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additional 2.2° alpha. The resulting 32.2° a would correspond to the corrected a at FS 13.56 
including wing induced angle of attack effects. The results are shown in figure 41 compared with 
the results presented above for a = 30° and /? = 0°. Although the results at the first station (where 
the estimated effect was only 1°) start to show good agreement, the change in the solution at 
FS 13.56 is not enough to attribute the differences to the neglected wing flowfield. Also note that if 
viscous effects had been included an even greater angle of attack change would be required to 
obtain agreement with experiment. Thus the wing flowfield effects do not entirely explain the poor 
agreement between experiment and computation at the suction peak locations. 

The flowfield results shown in figures 42 to 53 have been generated as in the case of F-5A 
using computational flow visualization techniques (PLOT3D). Figure 42 shows the inviscid 
surface oil flow pattern for a = 30° and p = 0°. The dye is injected at two different axial stations 
to observe the formation of separation lines, if any. The inviscid separation can be clearly seen to 
be occurring at the sharp edges as expected. The turbulent surface flow pattern is completely 
different, as shown in figure 43 for the same case. The primary separation lines start from the edge 
as in the case of inviscid flow, and an attachment line occurs at the centerline. However a 
secondary separation line is shown in the figure. 

Figures 44 and 45 show the cross sectional pressure contours at the axial station x = 13.56 
in. from nose for the inviscid and viscous flows respectively at a = 30° and (5 = 0°. The 
pressures are seen to be symmetrical about the plane of symmetry for both the inviscid and 
turbulent cases. The vortex pattern and symmetry are more clearly seen in the figures 46 and 47 
which are shaded based on the stagnation pressure. The turbulent flow also shows the secondary 
vortices being shed from under the primary vortices further downstream from the nose. 

Figures 48 to 53 show the similar flow visualization pictures for the case of a = 30° and 
(5 - 10°. Figure 48 shows the inviscid surface oil flow pattern. The asymmetry pattern is readily 
apparent. The lines that end abruptly arise when the streamline reaches the edge of the 
computational grid. After this photo session we learned how to operate PLOT3D to avoid this 
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premature termination of the oil flow line. However, the figure demonstrates that the attachment 
line has shifted due to sideslip. Figure 49 presents the results when viscous effects are included. 
Qualitatively the results are the same as in the zero sideslip case, but show the asymmetry 
introduced by the sideslip. Figures 50 and 51 show the pressure contours with the magnitudes of 
the pressures shown on the left. In both the turbulent and inviscid cases the leeward vortex is 
farther away from the surface than the windward vortex. Figures 52 and 53 show the stagnation 
pressure contours, and define the vortex location more precisely. 

A MODEL FOR FOREBODY DESIGN: THE GENERIC FOREBODY 


Passive tailoring of the aerodynamic properties of forebodies can only be achieved through 
the geometry of the forebody. As a means of beginning the study of geometry effects on forebody 
characteristics, a class of shapes possessing the ability to produce shapes of interest has been 
defined. In particular, the body can be considered to be composed of essentially independent 
cross-section and planform lines. This forebody makes use of the equation of a super-ellipse to 
obtain both cross section and planform lines. The super-ellipse can: 


• Recover a circular cross-section 

• Produce elliptical cross-sections 

• Capture the key characteristics of chine forebodies, such as the one tested by 
Erickson and being used on the YF-23. 

• Produce a straight sidewall at the max half breadth line with various slopes 

• Be generated by a simple analytical equation 
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The equation for the forebody cross section is: 




= 1 


where n and m adjustable coefficients. The constants a and b correspond to the max 
halfbreadth and centerlines respectively. Depending on the value of n and m , the equation can be 
made to meet all the requirements specified above. The case n = m = 0 corresponds to the 
standard ellipse, where the body is circular when a = b. When n = -1 the sidewall is linear at 
the max half breadth line, forming a distinct crease line. When n < -1 the body cross-section takes 
on the cusped or chine-like shape. The derivative oiztb with respect to y / a is: 



where z = z/b and y = y / a. As y 1, the slope becomes: 


dz 

dy 


<0 

-(2 + rn)y 1+m 


n>— 1 
n <- 1 

n = - 1 


Extreme flexibility is provided by allowing n, m, a, b to be functions of the axial distance x, 
and providing for different cross-sections to be used above and below the max half-breadth line. 
Notice that when n = -1 the value of m can be used to control the slope of the sidewall at the 
crease line. 

An example of the range of shapes available within this simple parametric geometry is 
illustrated in the figures. Figure 54 provides an upper quadrant section to illustrate the manner in 
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which the cross-section can change from a chine, through the straight line sidewall at the max 
halfbreadth to the elliptic and even "fuller" elliptic shape originally envisioned for super ellipses 
simply by changing the value of the exponent coefficient n. As n becomes very large the shape 
will approach a rectangular cross section. Figure 55 shows the special case where the wall slope is 
finite and non-zero at the max half-breadth (n = -1). In this case the value of the slope is 
controlled by the value of m. The figure demonstrates the range of control that is available, and 
the effect of sidewall slope selection on the entire cross-section. 

The entire cross section is shown for a series of axial station in figure 56. Figure 56a shows 
a standard elliptical cross section, with different upper and lower surface centerlines, while figure 
56b provides similar results for the case with the straight line side at the max half breadth, and 
figure 56c provides the equivalent view for the chine shaped case. The range of cross sectional 
shapes, together with the possibility of allowing them to vary character with axial distance provides 
an extremely broad design space to investigate aerodynamic tailoring of forebody characteristics 
through geometric design. 


CONCLUSIONS 

In this study the use of a Navier-Stokes code to obtain credible results has been demonstrated 
for high angle of attack forebody aerodynamics at combined angle of attack and sideslip.The F-5A 
forebody experimental results have been simulated computationally. The computed results nearly 
duplicated the experimental results. By comparing inviscid and viscous computational solutions the 
role of viscosity in creating the stabilizing effect of the forebody has been explicitly identified. 

When considering extreme forebody shapes such as the Erickson forebody, the code has also 
been shown to be capable of obtaining solutions in qualitative agreement with measured pressure 
distributions, although differences between computed and measured pressure distributions remain. 
Further analysis must be carried out to determine if grid resolution or turbulence model deficiencies 
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are responsible for the lack of quantitative agreement. 

In obtaining the solutions used to determine the directional stability derivatives it was found 
that the solution had to be “fully converged”. Because the derivative was obtained as the difference 
between two small numbers, the value of C n ^ only started to converge after around 4000 
iterations. 

An analytic generic forebody model that is extremely versatile yet very simple has been 
proposed to study the relation between forebody geometry and aerodynamic characteristics. 
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APPENDIX A 


MODIFICATIONS TO cfl3d ( version 1.1) 

The following modifications and additions were made to cf!3d ( version 1.1 ) 

Where : bc.f ; subroutine be 

Purpose : Periodic boundary condition at the grid interface of a wraparound grid used for 

flows with sideslip 

New : 

if (mtypej (1) .eq.2001 .or. mtype j (2) .eq.2001) 

. call bc2001 ( jdim, kdim, idim,w(lq) ,w(lqjO) ,w(lqkO) ,w(lqiO) , 

. w(lsj) ,w(lsk) ,w(lsi) ) 


Where : 
Purpose : 

New : 


bc.f ; completely new subroutine 

Periodic boundary condition at the grid interface of a wraparound grid used for 
flows with sideslip 

subroutine bc2001 (jdim, kdim, idim, q, qjO, qkO,qiO, s j, sk, si) 


common / reyue / i vise (3) ,reue,tinf,isnd, c2spe 
c 

common /info/ title(20) ,xmach, alpha, beta, nit, ntt,dt,fmax, idiag(3) , 
. nitfo, if lagts, rkap (3) , istrag, impl, iru, n ju, i jac, iaf , nres, iaf a, 

. levelb (5) ,mgflag, iconsf ,mseq, ncycl (5) , levelt (5) ,nitfol (5) ,mmx, 

. imesh,ngam,nsm(5) ,iflim(3) 
common/fluid/gamma, gml, gpl, gmlg, gplg, ggml 
common/ivals /pO , rhoO , cO , uO , vO , wO , et 0 , hO , pt 0 , qi v ( 5 ) 
common /te/ jtel, jte2,ktel,kte2,itel,ite2 
c 

common/bv/ibci (2) ,ibcj (2) ,ibck(2) ,ibcmi (5,2) ,ibcmj (5,2) ,ibcmk(5,2) 
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common /type/ mtype j (2) ,mtypek (2) ,mtypei (2) 
c 

common/unst/iunst , time, rf req, alf , alphau, cloc, rf reqO 
common /ski ton/ i ski ton 
c 

dimension q( jdim, kdim, idim, 5) ,qiO ( jdim,kdim, 5,4), 

qjO (kdim, idim-1, 5,4), qkO ( jdim, idim-1, 5,4) 
dimension sk ( jdim, kdim, idim-1, 5) , 

. si (jdim, kdim, idim, 5) , s j (jdim, kdim, idim-1, 5) 
dimension tsym(5) 

data tsym / I.e0,l.e0,-l.e0,l.e0,l.e0 / 
c 

jdiml=jdim-l 

kdiml=kdim-l 

idiml=idim-l 

c 

c 

c 

Q ★★★★★★★★★★★★★★★★★★★★ ★★★ ★ 

c * yawed flow boundary condition mtype=2001 * 

Q ******************** *** * 

c 

c j=l boundary flow across x-z plane mtype 2001 

Q* ************************************************* ************************ 

C 

c reflection - symmetry plane j=0/j=jdim 

c 

if (mtype j (1) .eq. 2001) then 
if (isklton.eq.l) write(ll,*) 

j=l boundary flow across x-z plane from j = jdim mtype 2001' 
do 38 i=l , idiml 
do 38 1=1,5 
do 38 k=l , kdiml 

qjO (k, i, 1, 1) = q( jdim-l,k, i, 1) 

qjO (k, i, 1, 2) = q(jdim-2,k,i,l) 

38 continue 
end if 
c 
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I********************************************************* ********** ******* 


c 

c j=jdim boundary flow across x-z plane mtype 2001 

q * ************* ************************************************** *★**★**★★* 

C 

c reflection - symmetry plane j=0/j=jdim 

c 

if (mtype j (2) .eq. 2001) then 
if (isklton.eq.l) write (11,*) 

. ' j=jdim boundary flow across x-z plane from j = 1 mtype 2001' 
do 39 i=l,idiml 
do 39 1=1,5 
do 39 k=l , kdiml 
qjO (k, i, 1, 3) = q(l,k,i,l) 

qjO (k, i, 1, 4) = q(2,k, l, 1) 

39 continue 
end if 
return 
end 

Where : cbsem.f ; subroutine plot3d 

Purpose : Calculate and print coefficient of pressure Cp instead of 

Old: 

write (17,29)i,j,k, x(j,k, i),y(j,k, i),z(j,k, i), 
q2,q3,q4,q5, tl, xml, pitot, edvis 

New : 

cp = 2.0/gamma/xmach**2* (q5 - 1.0) 

write (17, 29) i, j,k,x(j,k, i) ,y ( j, k, i) , z ( j, k, i) , 

. q2,q3,q4,cp,tl, xml, pitot, edvis 

Where : cbsem.f ; subroutine plot3d 

Purpose : At a comer point, for wraparound grids, this is more appropriate 

Old: 

if (k.eq.kdim) kd = kdiml 
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c 


if (j.eq.l .or. j.eq.jdim) then 
if (j.eq.jdim) jd = jdiml 
corner points 

xw( jw,kw,iw,l) = 0.5* (q(jd,kd,id, l)+q(jd,kd,idl,l) ) 


New : 

m=2 

if (k .eq. kdim ) m=4 
if (j .eq. 1) 

xw( jw,kw,iw,l) = .25* (qkO ( jd,id,l,m) +qk0 (jdiml, id, l,m) 

• + qkO ( jd, idl, l,m) + qkO (jdiml, idl, l,m) ) 

if (j .eq. jdim) 

xw ( jw, kw, iw, 1) = . 25* (qkO ( jd, id, l,m) +qk0 (1, id, l,m) 

+ qkO (jd, idl,l,m) + qkO(l,idl,l,m) ) 


Where : lbcx.f ; subroutine force 

Purpose : Calculate yawing moment 

Old: 


subroutine force (jdim, kdim, idim,x,y, z, sk, q, 

. cl, cd, cz, cm, chd, swet, iOO, ub, vb, wb, vmu, vol) 
c 

Q** ********************************************************************* 


c Purpose: Integrate the forces on the body. 


t -k-k-k-k-k 


c implicit half precision (a-h, o-z) 

common /fsum/ sref,cref,bref,xmc,ymc, zmc 
common/ ivals/pO, rhoO, cO, uO, vO, wO, etO, h0,pt0, qiv (5) 
dimension ub ( jdim*kdim, 1 ) , vb ( jdim* kdim, 1 ) , 

+wb ( jdim* kdim, 1 ) , vmu ( jdim-1 , 1 ) , vol ( jdim* kdim, 1 ) 
common /te/ jtel, jte2,ktel,kte2,itel,ite2 

common /info/ title (20) ,xmach, alpha, beta, nit, ntt,dt,fmax, idiag (3) , 
. nitfo, iflagts, rkap (3) , istrag, impl, iru, n ju, i jac, iaf , nres, iafa, 

. levelb (5) ,mgflag,iconsf,mseq,ncycl (5) ,levelt (5) ,nitfol (5) ,mmx. 
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. imesh,ngam, nsm(5) ,iflim(3) 

dimension x ( jdim, kdim, idim) , y ( jdim, kdim, idim) , z ( jdim, kdim, idim) 
dimension sk( jdim*kdim, idim-1, 5) ,q( jdim, kdim, idim, 5) 
common/fluid/ gamma, gml, gpl, gmlg, gplg, ggml 
common / reyue / ivisc(3) ,reue,tinf, isnd, c2spe 
c 

common/unst/iunst, time, rfreq, alf , alphau, cloc, rfreqO 
c 

al=alpha+alf 

c 

cpc=2 . eO/ (gamma*xmach*xmach) 
c 

cosa=cos (al) 
sina=sin(al) 
c 

cl=0.e0 
cd=0.e0 
cz=0.e0 
cm=0.e0 
chd= O.eO 
swet=0 . eO 
ist=itel 
ifn=ite2-l 
jte21=jte2-l 
const=4 . / (reue*xmach) 
do 9000 i=ist,ifn 
cxl=0.e0 
cyl=0.e0 
czl=0.e0 
cml=0.e0 
chdl=0.e0 
xas=0.e0 
yas=0.e0 
zas=0.e0 
swetl=0.e0 
jte21=jte2-l 
do 50 j=jtel,jte21 
c 

xa=.25e0*( x(j,l,i) + x(j+l,l,i) + x(j,l,i+l) + x(j+l,l,i+l) ) 
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ya=.25eO* ( y(j,l,i) + y(j+l,l,i) + y(j,l,i+l) + y(j+l,l,i+l) ) 

za=.25e0* ( z(j,l,i) + z(j+l,l,i) + z(j,l,i+l) + z(j+l,l,i+l) ) 

xas=xas+xa 

yas=yas+ya 

zas=zas+za 

dcp=- (q( j, 1, i, 5) /pO-l.eO) *cpc*sk ( j, i, 4) 
dcx=dcp* sk ( j , i , 1 ) 
dcy=dcp* sk ( j , i , 2 ) 
dcz=dcp*sk ( j,i, 3) 

if (ivisc(3) .gt.O) then 

tau=vmu ( j , i ) *const/vol ( j , i ) *sk ( j , i, 4 ) **2 

vnorm=ub (j,i)*sk(j,i,l) + vb (j,i)*sk(j,i,2)+wb(j,i)*sk(j,i,3) 

dcx=dcx+tau* (iab( j, i) -vnonn*sk ( j,i,l) ) 

dcz=dcz+tau* (wb(j,i)-vnorm*sk(j,i,3) ) 

dcy=dcy+tau* ( vb ( j , i ) -vnorm* sk ( j , i , 2 ) ) 

end if 

chdl=chdl+abs (sk( j, i, 3) ) *sk ( j, i, 4) 

swetl=swetl+sk( j,i, 4) 

cxl=cxl+dcx 

cyl=cyl+dcy 

czl=czl+dcz 

50 cml=cml-dcz* (xa-xmc) +dcx* (za-zmc) 
xas=xas/float ( jte2-jtel) 
yas=yas/float ( jte2-jtel) 
zas=zas/float ( jte2-jtel) 
cds= cxl*cosa+czl*sina 
cls=-cxl*sina+czl*cosa 
cms=cml 
chds=chdl 
swets=swetl 
cl=cl+cls 
cd=cd+cds 
cz=cz+czl 
cm=cm+ans 
chd=chd+chds 
swet=swet+swets 

if(ntt.le.l .or. iOO.eq.l) then 
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c 


if ( chds . le . 0 . ) chds=l . 
cls=2 . 0*cls/chds 
cds=2 . 0*cds/chds 
czl=2 . 0*czl/chds 
cyl=2 . 0*cyl/chds 
c 

if (i.eq.ist) write (11, 8613) 

8613 format (/30x,41hSUMMARY OF INCONSTANT K=1 SECTIONAL LOADS,/, 

. 23x, 54hSref-s is 1/2 projected area in X-Y plane, if positive) 
if (i.eq.ist) write (11, 8713) 

8713 format (/lx, 2x,lhl, 5x, 6hXavg-s, 6x, 6hYavg-s, 6x, 6hZavg-s, 6x, 

. 6hSref-s, 6x, 6hSwet-s, 7x, 4hCl-s, 8x, 4hCd-s, 8x, 4hCz-s, 

. 8x, 4hCy-s, 7x, 7hRe-Xavg) 
c 

xrep = reue*xas 

write (11, 1318) i,xas,yas, zas, 

. chds, swets, els, cds, czl, cyl,xrep 
1318 format (Ix,i3,10el2.3) 
c 

end if 

9000 continue 
return 
end 


New : 


subroutine force (jdim, kdim, idim,x,y, z, sk, q, 

. cl, cd, cz, cm, cx, cy, cn, chd, swet, iOO, ub, vb, wb, vmu, vol) 


c Purpose : Integrate the forces on the body . 

C*********************************************************************** 


c implicit half precision (a-h,o-z) 

common /fsum/ sref, cref,bref,xmc,ymc, zme 
Common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0,qiv(5) 
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c 


c 


c 


c 


c 


dimension ub ( jdim*kdim, 1 ) , vb ( jdim* kdim, 1 ) , 

+wb ( jdim* kdim, 1) , vmu ( jdim-1 , 1 ) , vol ( jdim* kdim, 1) 
common /te/ jtel, jte2,ktel,kte2,itel,ite2 

common /info/ title(20) ,xmach, alpha, beta, nit, ntt,dt,fmax, idiag(3) , 
. nitfo, if lagts, rkap (3) , istrag, impl, iru, n ju, i jac, iaf , nres, iaf a, 

. levelb(5) ,mgflag, iconsf,mseq, ncycl (5) , levelt (5) , nitfol (5) ,mmx, 

. imesh,ngam, nsm(5) ,iflim(3) 

dimension x ( jdim, kdim, idim) ,y ( jdim, kdim, idim) , z ( jdim,kdim, idim) 
dimension sk ( jdim* kdim, idim-1, 5) , q( jdim, kdim, idim, 5) 
common/fluid/gamma, gml, gpl, gmlg, gplg, ggml 
common / reyue / ivisc (3) , reue, tinf , isnd, c2spe 

common/unst/iunst, time, rfreq, alf , alphau, cloc, rfreqO 

al=alpha+alf 

cpc=2 . eO/ (gamma*xmach*xmach) 

cosa=cos (al) 
sina=sin(al) 

cl=0.e0 

cd=0.e0 

cx=0.e0 

cy=0 . eO 

cz=0 .eO 

cm=0.e0 

cn=0.e0 

chd= O.eO 

swet=0 . eO 

ist=itel 

ifn=ite2-l 

jte21=jte2-l 

const=4 . / (reue*xmach) 

do 9000 i=ist, ifn 

cxl=0.e0 

cyl=0 . eO 

czl=0.e0 

cml=0.e0 

cnl=0.e0 
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chdl=0.e0 

xas=0.e0 

yas=0.e0 

zas=0.e0 

swetl=0.e0 

jte21=jte2-l 

do 50 j=jtel,jte21 

xa=.25e0*( x(j,l,i) + x(j+l,l,i) + x(j,l,i+l) + x(j+l,l,i+l) ) 

ya=. 25e0* ( y(j,l,i) + y(j+l,l,i) + y(j,l,i+l) + y(j+l,l,i+l) ) 

za=. 25e0* ( z(j,l,i) + z(j+l,l,i) + z(j,l,i+l) + z(j+l,l,i+l) ) 

xas=xas+xa 

yas=yas+ya 

zas=zas+za 

dcp=- (q( j, 1, i, 5) /p0-l.e0) *cpc*sk ( j,i, 4) 
dcx=dcp* sk ( j , i , 1 ) 
dcy=dcp*sk ( j, i, 2) 
dcz=dcp*sk ( j, i, 3) 

if (ivisc (3) .gt.O) then 
tau=vmu(j,i)*const/vol( j,i)*sk(j,i,4)**2 

vnorm=ub(j,i)*sk(j,i,l)+vb(j,i)*sk(j,i,2)+wb(j,i)*sk(j,i,3) 

dcx=dcx+t au* (ub ( j, i) -vnorm* sk ( j , i , 1 ) ) 

dcz=dcz+tau* (wb ( j, i) -vnorm*sk ( j, i, 3) ) 

dcy=dcy+tau* ( vb ( j , i ) -vnorm* sk ( j , i , 2 ) ) 

end if 

chdl=chdl+abs(sk(j,i,3) )*sk(j,i,4) 

swet l=swet 1+ sk ( j , i , 4 ) 

cxl=cxl+dcx 

cyl=cyl+dcy 

czl=czl+dcz 

cnl=cnl-dcy* (xa-xmc) +dcx* (ya-ymc) 

50 cml=cml-dcz* (xa-xmc) +dcx* (za-zmc) 
xas=xas/float ( jte2-jtel) 
yas=yas/float ( jte2- jtel) 
zas=zas/float ( jte2- jtel) 
cds= cxl*cosa+czl*sina 
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cls=-cxl*sina+czl*cosa 

cms=cml 

chds=chdl 

swets=swetl 

cl=cl+cls 

cd=cd+cds 

cx=cx+cxl 

cy=cy+cyl 

cz=cz+czl 

cm=cm+cms 

cn=cn+cnl 

chd=chd+chds 

swet=swet+ swet s 

if(ntt.le.l .or. iOO.eq.l) then 
c 

if ( chds . le . 0 . ) chds=l . 
cls=2 . 0*cls/chds 
cds=2 . 0*cds/chds 
czl=2 . 0*czl/chds 
cyl=2 . 0*cyl/chds 
c 

if(i.eq.ist) write (11, 8613) 

8613 format (/30x,41hSUMMARY OF I=CONSTANT K=1 SECTIONAL LOADS,/, 

. 23x, 54hSref-s is 1/2 projected area in X-Y plane, if positive) 
if(i.eq.ist) write (11, 8713) 

8713 format (/lx,2x,lhl, 4x, 6hXavg-s, 5x, 6hYavg-s, 5x, 6hZavg-s, 5x, ' 

. 6hSref-s, 5x, 6hSwet-s, 6x, 4hCl-s, 7x, 4hCn-s, 7x, 4hCx-s, 7x, 4hCz-s, 

. 7x, 4hCy-s, 6x, 7hRe-Xavg) 
c 

xrep = reue*xas 

write (11, 1318) i,xas,yas, zas, 

. chds, swets, els, cnl, cxl, czl, cyl, xrep 
1318 format (lx, i3, Hell. 3) 
c 

end if 

9000 continue 
return 
end 
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Where : 
Purpose 

Old: 


c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


rhs.f ; subroutine resp 

Calculate and print block summary of axial, normal and side forces and also of 
yawing moment 


call force ( jdim, kdim, idim, x,y, z, sk,q, 

. cl, cd, czz, cm, chd, swet, icall,q(l, 2) ,q(l,3),q(l,4) ,vmu,vol) 


clt 

= cl/sref 

cmt 

= cm/ (sref*cref ) 

cdt 

= cd/sref 


print block summary 

if (icyc.le.l .or. icyc.eq.ncyc) then 
write (11, 677) nbl 

677 format (/10x, 37hSUMMARY OF FORCES AND MOMENTS - BLOCK, i5) 
write (11, 603) 

603 format (/lx, llhwetted area, 2x, 14hX-Y proj. area, 

5x, 4hCL-b, lOx, 4hCD-b, lOx, 4hCM-b) 
write (11, 653) swet, chd, clt, cdt, cmt 
653 format (lx,el2.5, 4 (2x,el2. 5) ) 
end if 

s» contributions of blocks on global level in clw, cdw 

if (level. eq.lglobal) then 
chdgp = chdgp + chd 

swetgp = swetgp + swet 

clw(nres) = clw(nres) + clt 

cmgp = cmgp + cmt 

cdw(nres) = cdw(nres) + cdt 

end if 

print summary of global blocks 

if (icyc.le.l .or. icyc.eq.ncyc) then 
if (nbl .eq. lblock) then 
write (11, 777) 
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777 format (/10x, 4 9hSUMMARY OF FORCES AND MOMENTS - ALL GLOBAL BLOCKS) 
write (11, 703) 

703 format (/lx, llhwetted area,2x, 14hX-Y proj. area, 

6x, 2hCL, 12x, 2hCD, 12x, 2hCM) 

write (11, 653) swetgp,chdgp,clw(nres) ,cdw(nres),cmgp 
end if 
end if 
670 continue 
c 

irite=l 

if (irite .eq. 1 .and. iunst .gt. 0) then 
call prntcp ( jdim,kdim, idim,wk,q) 
end if 
c 

c print residual and lift 

c 

if (icyc.eq.l .or. (icyc.eq.2 .and. nbl . eq. lblock) 

. .or. icyc.eq.ncyc) then 
if (iunst .It. 2) then 
write (11, 5002) 

5002 format (/lx, 5hlevel, lx, 5hblock, lx, 9hiteration, 3x, 

. 8hresidual,4x,10htotal res.,3x,7hrkap(i) ,3x, 7hrkap(j) , 

. 3x, 7hrkap (k) , 3x,10hlift coef . , 3x, lOhdrag coef.) 

c 

else 

write (11, 4682) 

4682 format (/lx, 5hlevel, lx, 5hblock, lx, 9hiteration, 3x, 

8hresidual,4x,10htotal res. , 3x, 7hrkap (i) , 3x, 7hrkap ( j) , 

3x, 7hrkap (k) , 3x, lOhlift coef ., 3x, lOhdrag coef., 

4x, 4htime, 6x, 5halpha) 

end if 
end if 

if (icyc.eq.l .or. (icyc.ge.2 .and. nbl. eq. lblock) 

. .or. icyc.eq.ncyc) then 
if (iunst .It. 2) then 

write (11, 5001) level, nbl, nres,rmst,rms (nres) ,rkap, 

. clw (nres) , cdw (nres) 

5001 format (lx, i5, lx, i5, 3x, i5, 2x, el2 . 5, lx, el2 . 5, 2x, f 6 . 2, 

. 4x,f6.2,4x,f6.2,3x,el2.5,lx,el2.5) 
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else 

alot= (alpha+alf ) *57 .2958 

write (11, 4683) level, nbl,nres,rmst,rms (nres) ,rkap, 

. clw (nres) ,cdw (nres) , time, alot 

4683 format (lx, i5,lx, i5, 3x,i5, 2x,el2. 5, lx,el2. 5, 2x,f6.2, 

. 4x,f6.2,4x,f6.2,3x,el2.5,lx,el2.5,lx,f9.4,lx,ell.4) 

end if 
end if 
c 

if (nbl.eq.lblock) then 


chdgp 

= 0 . 

swetgp 

= 0 . 

cmgp 

= 0 . 

end if 


67 continue 



New : 

call force ( jdim, kdim, idim,x,y, z, sk,q, cl, cd, cz, cm, 

. cx, cy, cn, chd, swet,icall,q(l,2) ,q(l,3) ,q(l,4) ,vmu,vol) 
c 


clt 

= cl/sref 

emt 

= cm/ (sref*cref ) 

cdt 

= cd/sref 

ext 

= cx/sref 

czt 

= cz/sref 

cyt 

= cy/sref 

ent 

= cn/ (sref*bref ) 


c 

c print block summary 

c 

if (icyc.le.l .or. icyc. eq. ncyc) then 
write (11, 677) nbl 

677 format (/10x, 37hSUMMARY OF FORCES AND MOMENTS - BLOCK, i5) 
write (11, 603) 

603 format (/lx, llhwetted area, 2x, 14hX-Y proj. area, 4x, 4hCL-b, 

. 9x, 4hCD-b, 9x, 4hCM-b, 9x, 4hCX-b, 9x, 4hCY-b, 9x, 4hCZ-b, 9x, 4hCN-b) 
write (11, 653) swet, chd, clt, cdt, ant, ext, cyt, czt, ent 
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653 format (lx,el2.5, 8 (2x,ell. 4) ) 
end if 
c 

c sum contributions of blocks on global level in clw, cdw 

c 

if (level. eq.lglobal) then 


chdgp 

= 

chdgp 

+ 

chd 

swetgp 

= 

swetgp 

+ 

swet 

clw (nres) 

= 

clw (nres) 

+ 

clt 

cmgp 

= 

cmgp 

+ 

cmt 

cdw (nres) 

= 

cdw (nres) 

+ 

cdt 

cxgp 

= 

cxgp 

+ 

ext 

czgp 

= 

czgp 

+ 

czt 

cygp 

= 

cygp 

+ 

cyt 

cngp 

= 

cngp 

+ 

ent 


end if 
c 

c print summary of global blocks 
c 

if (icyc.le.l .or. icyc.eq.ncyc) then 
if(nbl .eq. lblock) then 
write (11, 777) 

777 format (/10x, 49hSUMMARY OF FORCES AND MOMENTS - ALL GLOBAL BLOCKS) 
write (11,703) 

703 format (/lx, llhwetted area, 2x, 14hX-Y proj. area, 

5x, 2hCL, llx, 2hCD, llx, 2hCM, llx, 2hCX, llx, 2hCY, 
llx, 2hCZ, llx, 2hCn) 

write (11, 653) swetgp, chdgp, clw (nres) , cdw (nres) , cmgp, cxgp, cygp, 

. czgp, cngp 

end if 
end if 
670 continue 
c 

irite=l 

if(irite .eq. 1 .and. iunst .gt. 0) then 
call prntcp( jdim,kdim,idim,wk,q) 
end if 
c 

c print residual and lift 
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if (icyc.eq.l .or. (icyc.eq.2 .and. nbl . eq. lblock) 

. .or. icyc.eq.ncyc) then 
if(iunst .It. 2) then 
write (11, 5002) 

5002 format (/lx, 5hlevel, lx, 5hblock, lx, 9hiteration, 3x, 

. 8hresidual, 4x, lOhtotal res. , 3x, 7hrkap (i) , 3x, 7hrkap ( j) , 

. 3x, 7hrkap(k) , 3x,10hlift coef . , 3x, lOhdrag coef., 

. 3x,7h cy ,3x, 7hcn ) 

c 

else 

write (11, 4682) 

4682 format (/lx, 5hlevel, lx, 5hblock, lx, 9hiteration, 3x, 

. 8hresidual, 4x, lOhtotal res. , 3x, 7hrkap (i) , 3x, 7hrkap ( j) , 

. 3x, 7hrkap (k) , 3x, lOhlift coef ., 3x, lOhdrag coef., 

. 4x, 4htime, 6x, 5halpha) 

end if 
end if 

if (icyc.eq.l .or. (icyc.ge.2 .and. nbl. eq. lblock) 

. .or. icyc.eq.ncyc) then 
if(iunst .It. 2) then 

write (11, 5001) level, nbl, nres, rmst,rms (nres) ,rkap, 

. clw (nres) ,cdw (nres) ,cygp,cngp 

5001 format (lx, i5,lx,i5, 3x,i5, 2x, el2. 5, lx, el2 . 5, 2x, f 6. 2, 

4x, f 6 . 2, 4x, f 6. 2, 3x, el2 . 5, lx, el2 . 5, lx, el2 . 5, lx, el2 . 5) 

else 

alot= (alpha+alf ) *57 .2958 

write (11, 4683) level, nbl, nres, rmst, rms (nres) ,rkap, 

. clw (nres) , cdw (nres) , time, alot 

4683 format (lx, i5, lx, i5, 3 x, i5, 2x, el2 . 5, lx, el2 . 5, 2x, f 6. 2, 

4x, f 6 . 2, 4x, f 6 . 2, 3x, el2 . 5, lx, el2 . 5, lx, f 9 . 4, lx, ell . 4 ) 

end if 
end if 
c 

if (nbl. eq. lblock) then 
chdgp = 0 . 

swetgp = 0 . 

cmgp = 0. 

cxgp = 0 . 

cygp = 0 . 
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Where : 
Purpose : 

Old: 

New : 

Where : 
Purpose : 

Old: 

New : 


czgp = 0. 

cngp = 0 . 

end if 
67 continue 


rhs.f ; subroutine turb 

Correction to the value of parameter cwk used in Baldwin Lomax model ( AIAA 
78-257 ) which was incorrect 

data aplus/26. e0/, ccp/1 . 6e0/, ckleb/ . 3e0/, cwk/ . 25e0/, vk/ . 4e0/, 
clauser/ . 0180e0/ 


data aplus/26 .e0/ , ccp/1 . 6e0/, ckleb/ . 3e0/, cwk/1 . OeO/, vk/. 4e0/, 
. clauser/. 0180e0/ 

rhs.f ; subroutine turb 

Change initialized value of fblmax to avoid division by zero. More of a 
precautionary measure as suggested by Dr. R.W.Walters 


fblmax = 0. 


fblmax = l.e-10 
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APPENDIX B 


F-5A FOREBODY GEOMETRY MATH MODEL 


The F-5 A forebody math model was constructed using conic lofting techniques ( ref 7 ) based on 
lofting information provided by Tom Heglund of Northrop, and verified by comparison with an 
actual nose provided by Sue Grafton. The computational model was extended to 225 in. (without 
the canopy) which was beyond the length of the wind tunnel model (175 in. ). 


program fSagrid ( F-5A forebody math model ) 
c F-5A forebody surface geometry and grid definition. 

c R.Ravi f August 1989 

c program generates the control and maximum half breadth points at 

c various cross sections as suggested by Tom Heglund of Northrop 

c and then generates the cross sectional shape using conic lofting 

c techniques (ref. 24 ) 

c 

c It also generates the F-5A surface grid and uses this 

c information to generate the input to the grid generator 

c which in our case was GRIDTOOL 

c 

c The length of the computational model was 225 in. (canopy 

c not included) 

c 

c It scales down the surface co-ordinates to compare with 

c Sue Graftons F-5A wind tunnel model 

c 

c It finally writes down the surface co-ordinates in PL0T3D format 

c 

c 

dimension x (100,10) ,y (100) , z (100, 10) ,xl (100,100) ,x2 (100,100) , 

1 zl (100,100) ,z2 (100,100) , xxx (100, 100,1) ,yyy (100, 100, 1) 

2 , zzz (100, 100, 1) , xx2 (100,100) , zz2 (100, 100) ,yext (15) 
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3 , slopel (100, 100) ,curv(100, 100) , s (100) 

c open (7, f ile= ' f 5asurface.out ' ) 

open (10, f ile= ' f 5asurface .plot3d' ) 

open (11, f ile=' f 5asurf ace .compare ' ) 

open (13, file='test' ) 

open (14, f ile= ' slopes ' ) 

open (15, f ile= ' curv' ) 

open (12, f ile= ’ f body . gridgen ' ) 

open (16, file= ' report ' ) 

pi=acos (-1 . ) 

C 

C IEXT = NUMBER OF AXIAL STATIONS AHEAD OF THE NOSE 

C INCLUDING THE NOSE STATION 

C JDIM NUMBER OF POINTS FROM THE TOP CENTER LINE 

C TO THE POINT WHERE THE FLAT SURFACE STARTS 

C Y AHEAD = DISTANCE AHEAD OF THE NOSE FROM WHERE 

C SUBSONIC GRID STARTS 

C 

jdim=25 
jdiminit= jdim 
iext=10 
j2= jdim*2-l 
j2init= j2 
yahead=100 . 0 
y (1) =42 .50 
y (2) =42 . 55 
y (3) =42 . 60 
y (4) =42 . 65 
y (5) =42 .70 
y (6) =42 .75 
y (7) =42 . 80 
y (8) =42 . 85 
y (9) =42 . 90 
y (10) =42 . 95 
y (11) =43 .00 
y (12) =43 . 50 
y (13) =44 . 00 
y (14) =44 .50 
y (15) =45 .00 
y (16)=45.50 


44 



y (17) =46 . 00 
y (18) =47 .00 
y (19) =48 . 00 
y (20) =49 . 00 
y (21) =50 • 00 
y (22) =52 .00 
z (1, 1) =-8 .589 
z (1, 6) =-8 .589 
z (2, 1) =-8 . 464 
z (2, 6) =-8 .745 
z (3, 1) =-8 . 408 
z (3, 6) =-8 .811 
z (4, 1) =-8 .364 
z (4, 6) =-8 . 861 
z (5,1) =-8 . 326 
z (5, 6) =-8 . 904 
z (6, 1) =-8 .291 
z (6, 6) =-8 . 942 
z (7, 1) =-8 .259 
z (7, 6) =-8 . 977 
z (8, 1) =-8 .230 
z (8, 6) =-9 . 008 
z (9, 1) =-8 .201 
z (9, 6) =-9 . 038 
z (10, 1) =-8 . 175 
z (10, 6) =-9 . 066 
z (11, 1) =-8 .149 
z (11, 6) =-9 . 092 
z (12, 1) =-7 . 929 
z (12, 6) =-9 .308 
z (13, 1) =-7 .745 
z (13, 6) =-9 . 477 
z (14, 1) =-7 .581 
z (14, 6) =-9 . 621 
z (15, 1) =-7 . 428 
z (15, 6) =-9 .749 
z (16, 1) =-7 .283 
z (16, 6) =-9 . 8 67 
z (17, 1) =-7 . 145 
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2(17, 6)=-9. 976 
z (18, 1) =-6 . 884 
z (18, 6) =-10 . 175 
z (19, 1) =-6 . 636 
z (19, 6) =-10 . 356 
z (20, 1) =-6 . 4 
z (20, 6) =-10 . 523 
z (21, 1) =-6 . 171 
z (21, 6) =-10 . 68 
z (22, 1) =-5 . 721 
z (22, 6) =-10 . 98 
x (1, 3) =0 . 0 
x (2, 3) =-0 . 175 
x (3, 3) =- .251 
x (4, 3) =- . 311 
x (5, 3) =-.363 
x (6, 3) =- . 409 
x( 7, 3)=-. 452 
x (8, 3) =- . 491 
x ( 9, 3) =- . 528 
x (10 , 3) =- .564 
x (11, 3) =- . 598 
x (12, 3) =-0 . 884 
x ( 13, 3) =-l . 118 
x (14, 3) =-l . 32 6 
x (15, 3) =-l . 517 
x (16, 3) =-l .696 
x (17, 3) =-l .866 
x (18, 3) =-2 . 185 
x (19, 3) =-2 . 484 
x (20, 3) =-2 .768 
x (21, 3) =-3 . 041 
x (22, 3) =-3 . 577 
z (1, 3) =-8 . 589 
z (22, 3) =-8 .147 
do 3 i=2, 21 

slope= (z (22, 3) -z (1, 3) ) / (y (22) -y (1) ) 
z (i, 3) =z (1, 3) + (slope* (y (i) -y (1) ) ) 

3 continue 

do 5 1=1,22 
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x (i, 1) =0 . 0 
x (i, 6) =0 . 0 
continue 
idim=ll 

iidim=22+idim-l 

y (22) =52 . 0 

y (iidim) =125 .00 

dely= (y (iidim) -y (22) ) / (idim-1) 

do 10 i=22, iidim 

x (i, 1) =0 . 0 

x (i, 6) =0 . 0 

yl= .01* (y (i) -38 . 5) 

al=-. 0 88 15+.235*yl-. 04904* (yl**2. ) 

a2=. 0457594* (yl**3 .) -0 . 0294590* (yl**4 .)+. 0077218* (yl**5.0) 

a3=-0 . 000571* (yl**6.) 

z(i,l)=(al+a2+a3)/.01 

bl=-.08815-.17125*yl+. 068785* (yl**2 .)+. 1089634* (yl**3.) 
b2=-. 1823104* (yl**4.)+. 110555* (yl**5 .)-. 0242 660* (yl**6.) 
z (i, 6) = (bl+b2) / . 01 
y3= (y (i) -42 .5) *0.01 

a4=.0065+.213*y3-. 269652* (y3**3) +. 3669951* (y3**4) - .25149* (y3**5) 
a5=. 0874987* (y3**6 .)-. 0120677* (y3**7.) 
x (i, 2 ) =-l . 0* (a4+a5) / . 01 

if ( (y (i) .ge. 52.0) .and. (y(i) .le. 70.0)) then 
z (i, 2) =-15.551+. 1746575*y(i) 
x(i,5)=-1.0*(-7.37+.191*y(i) ) 
x (i, 3) =-l . 0* (-10. 35644+. 267 94 92*y(i) ) 
x (i, 4) =x (i, 3) 

else 

y4=.01*(y(i)-70.0) 

z (i, 2) =(-.03325+.1746575*y4+. 0215199* (y4**3. ))/. 01 
b4=.06+.191*y4-l. 1941 603* (y4**3 .) +3 . 5325676* (y4**4.) 
b5=-4. 7103427* (y4**5) +2.9985414* (y4**6) -. 7374221* (y4**7) 
x (i, 5) =-l . * (b4+b5) / . 01 

el=. 08 4+. 267 9492*y4-2. 0762353* (y4**3) +8 . 9749473* (y4**4) 
e2=-19. 2071906* (y4**5.) +20. 6116655* (y4**6.) 
e3=-8. 8061808* (y4**7.) 
x (i, 3) =-l .* (el+e2+e3) / . 01 
x(i f 4) =x (i, 3) 
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endif 

y5=. 01* (y (i) -42.5) 

b6=-. 09175- . 096*y5+. 3283568* (y5**3.)-. 8070199* (y5**4.) 
b7=. 8932135* (y5**5 4656938* (y5**6 .)+. 0929695* (y5**7 . ) 
z (i, 5) = (b6+b7) / . 01 
z (i, 3) =-11 . 0949+. 0589595*y (i) 
z (i, 4) =z (i, 3) 
y (i+1) =y (i) +dely 
10 continue 

idiml=7 

iidiml=iidim+idiml-l 
y (iidim) =125 . 00 
y (iidiml) =175 . 00 

dely= (y (iidiml) -y (iidim) ) / (idiml-1) 
do 20 i=iidim, iidiml 
x (i, 1) =0 . 0 
x (i, 6) =0 . 0 
yl=.01*(y(i)-125.) 

cl= . 0 9505 6+. 1965665*yl-. 017364* (yl**2 .)-. 1030717* (yl**3.) 
c2=. 8564063* (yl**4 .) -2 . 0571887* (yl**5 .) +1 . 5294378* (yl**6.) 
z (i, 1) = (cl+c2) / .01 
y2=.01*(y(i)-38.5) 

dl=- . 08815- . 17125*y2+ .068785* (y2**2 .)+. 1089634* (y2**3.) 
d2=-. 1823104* (y2**4.)+. 110555* (y2**5 .)-. 024266* (y2**6.) 
z (i, 6) = (dl+d2) / . 01 
y3= (y (i) -42 .5) *0.01 

b4=. 0065+.213*y3-.269652* (y3**3)+. 3669951* (y3**4) - .25149* (y3**5) 
b5=. 0874987* (y3**6 .)-. 0120677* (y3**7.) 
x(i,2)=-l.*(b4+b5)/0.01 
y4=.01*(y(i)-125.0) 

d3=. 066392+ . 1 9418 68*y4+. 0645597* (y4**2 .)-. 3888 642* (y4**3.) 
d4=2. 6147284* (y4**4.) -6. 9406616* (y4**5.) +5. 5973024* (y4**6.) 
z (i, 2) = (d3+d4) / . 01 
y5=.01*(y(i)-70.0) 

b6=.06+.191*y5-l. 1941603* (y5**3.) +3. 532567 6* (y5**4.) 
b7=-4. 7103427* (y5**5 .) +2 . 9985414* (y5**6 .)-. 7374221* (y5**7 . ) 
x (i, 5) =-l . * (b6+b7) / . 01 
y6= . 01* (y (i) -42 .5) 

b8=-. 09175- . 096*y6+. 3283568* (y6**3.)-. 8070199* (y6**4.) 
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b9=. 8 932135* (y6**5. )-. 4656938* <y6**6.)+. 0929695* (y6**7.) 
z(i,5) = (b8+b9)/.01 
y7=. 01* (y (i) -72 . 3) 

el=.09025+. 2 679492*y7-. 252125* (y7**2 .)+. 1235176* (y7**3.) 

e2=-. 0352414* (y7**4 .)+. 005155* (y7**5. ) 

x (i, 3) =-l . * (el+e2) / . 01 

x(i,4)=x(i,3) 

yl= . 01* (y (i) -125 . ) 

e3=-. 03725+ . 058 9595*yl+13. 7487102* (yl**3 . ) -98 . 7 915275* (yl**4.) 
e4=300. 6858408* (yl**5 .) -430 . 699982* (yl**6 .) +236 . 8200513* (yl**7 . ) 
z (i, 3) = (e3+e4) / . 01 
if (y (i) .le. 150.) then 
yl= .01* (y (i) -125 . ) 

e5=- . 03725+ . 058 9595 *yl- .05* (yl**2. ) -8.9138868* (yl**3.) 
e6=65. 1205618* (yl**4.) -191.1452871* (yl**5.) 
el =201 . 4735881* (yl**6. ) 
z (i, 4) = (e5+e6+e7) / . 01 

else 

y8=. 01* (y (i) -150 . ) 

e8=- . 04655-. 085*y8- .025* (y8**2 . ) +1 . 1163* (y8**3.) 
e9=-16 . 0196* (y8**4 . ) +79 . 9248* (y8**5 . ) -121 . 056* (y8**6.) 
z (i, 4) = (e8+e9) / . 01 

endif 

y ( i+1 ) =y ( i ) +dely 
20 continue 

idim2=6 

iidim2=iidiml+idim2 
y (iidiml) =175 . 00 
y(iidim2) =225.0 

dely=(y (iidim2) -y (iidiml) ) / (iidim2 -iidiml) 

do 700 i=iidiml, iidim2 

x(i, 1) =0 . 0 

yl0= . 01* (y (i) -175 . 0) 

gl=. 18 925+. 174*yl0-. 025* (yl0**2. ) -.0891784* (yl0**3.) 

g2=-. 404712* (yl0**4 .) +1 . 3978099* (yl0**5 .) -1 . 4735205* (yl0**6 . ) 

g3=. 5163509* (yl0**7.) 

z (i, 1) = (gl+g2+g3) / . 01 

x (i, 6) =0 . 0 

yll= (y (i) -175 . ) * . 01 

g4=- . 1826+ . 025* (yll**2.)+. 000825* (yll**3.) 
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g5=-. 00451* (yl 1**4.) _ . 00006* (yll**5.) 
z(i / 6)=(g4+g5)/.01 

if ( y(i) .ge. 175. .and. y(i) .le. 220.) then 
yl2=.01*(y(i)-42.5) 

g6=.0065+.213*yl2-. 2 696520* (yl2**3 .)+. 3669951* (yl2**4 . ) 
g7=-. 25149* (yl2**5)+. 0874987* (yl2**6)-. 0120677* (yl2**7) 
x (i, 2) =-l . * (g6+g7) / . 01 
yl3=. 01* (y (i) -72 .3) 

gl0=.09025+.267949*yl3-. 25212* (yl3**2) +. 1235176* (yl3**3) 
gll— . 0352414* (yl3**4.)+. 005155* (yl3**5.) 
x (i, 3) =-l . * (glO+gll) / .01 
x (i, 4) =x(i, 3) 

else 

x(i,2)=-15.50 
x (i, 3) =-20.25 
x (i, 4) =x (i, 3) 

endif 

g8= . 165+ . 155*yll- .025* (yll**2) -. 311936* (yll**3) +. 826826* (yll**4) 
g9=— . 9010034* (yll**5 .)+. 3592734* (yll**6 .)-. 01816* (yll**7.) 
z(i f 2)=(g8+g9)/.01 

gl2= . 05325+ . 085*yll- .075* (yll**2 ) + . 1473* (yll**3) -. 33543* (yll**4) 
gl3=. 5778333* (yll**5.) -.5472592* (yll**6 .)+. 1940741* (yll**7.) 
z (i, 3) = (gl2+gl3) / . 01 

if (y(i) .ge.175. .and. y(i) .le. 200.) then 

gl4=-.066-.0377*yll+6. 3274794* (yll**3) -54.34748* (yll**4) 
gl5=183. 0187488* (yll**5. ) -221.217856* (yll**6.) 
z (i, 4) = (gl4+gl5) / . 01 
gl6=.141+.016*yll-. 04* (yll**2) 

1 -. 6661333* (yll**3) +5.7792* (yll**6) 

gl7=-17 .80224* (yll**5 . ) +19 . 6949333* (yll**6 . ) 
x (i, 5) =-l . * (gl6+gl7) / . 01 

else 

z (i, 4) =-12. 4874+. 0303714 *y(i) 
x (i, 5) =-14 .209 

endif 

z (i, 5) =z (i, 6) +3 . 45 
y (i+1) =y (i) +dely 
700 continue 

do 100 i=l, 22 
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am=x ( i, 3 ) 
an=z (i, 3) 
au=x (i, 1) 
av=z (i, 1) 
ag=x(i, 6) 
ah=z (i f 6) 
aa=au-am 
ab=av-an 
ac=ag-am 
ad=ah-an 

c using a parabola for sections before y(i)<50.0 

pup=-l . 0* (ab/aa) 
sup= (4.0* (ab**2 . ) ) /aa 
pdown=-l . 0* (ad/ac) 
sdown= (4.0* (ad**2 . ) ) /ac 
do 105 j=l, jdiminit 

c xl (i, j) =aa*cos ( ( ( j-1) *pi) / (2 . 0* ( jdiminit-1) ) ) 

xl (i, j) =aa-aa*sin ( ( (j-1) *pi) / (2.0* (jdiminit-1) ) ) 
if (i .eq. 1) then 
zl (i, j ) =0 . 0 
go to 105 

endif 

zl (i, j) =pup*xl (i, j) + ( (abs (sup*xl (i, j) ) ) **0 . 5) 

105 continue 

do 110 j= jdiminit , j2init 
j j j =2* jdiminit- j 
xl (i f j ) =xl (if j j j) 
if (i .eq. 1) then 
zl (1, j)=0.0 
go to 110 

endif 

zl (if j) =pdown*xl (i f j) - ( (abs (sdown*xl (i, j) ) ) **0.5) 
110 continue 

do 115 j=lf j2 
x2 (if j) =xl (if j) +am 
z2 (i, j)=zl (if j) +an 
115 continue 

do 102 j=lf j2init 

write (13f 103) i f j f xl (i, j) f zl (if j) f -x2 (i, j) , z2 (i f j) 
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103 format (lx, 2 (i3, 2x) , lx, 4 (f 14 . 4, lx) ) 

102 continue 

100 continue 

do 200 i=23, 44 
if (i .ge. 33) then 
jdim=21 
j2= (2* jdim) -1 
jdimf inal= jdim 
j2f inal= j2 

endif 

if ((i .ge.23) .and. (i .le. 32)) then 
jdim=22 
j2= (2* jdim) -1 
j dimint er= jdim 
j2inter= j2 

endif 
aml=x (i, 3) 
anl=z (i, 3) 
am2=x (i, 4) 
an2=z (i, 4) 
aul=x (i, 1) 
avl=z (i, 1) 
au2=x (i, 2) 
av2=z (i, 2) 
agl=x(i, 6) 
ahl=z (i, 6) 
ag2=x (i, 5) 
ah2=z (i, 5) 
au3=aul-aml 
av3=avl-anl 
au4=au2-aml 
av4=av2-anl 

ak« ( (av4- (av3*au4/au3) ) **2 . ) / (au4* (av3-av4) ) 
as=ak*av3 

ar=ak* ( .25*ak- (av3/au3) ) 
ap=-0 .25*ak- ( ,25*ak- (av3/au3) ) 
do 205 j=l,jdim-l 

xl (i, j ) =au3-au3*sin ( ( ( j-1) *pi) / (2 . 0* < jdim-2) ) ) 
zl (i, j ) =ap*xl (i, j) + ( (abs (ar* (xl (i, j) **2 . ) +as*xl (i f j) ) ) **0.5) 
205 continue 


52 



do 210 j=l,jdim-l 
x2 ( i , j ) =xl ( i , j ) +aml 
z2 (i, j) =zl (i, j) +anl 
210 continue 

ag3=agl-am2 

ah3=ahl-an2 

ag4=ag2-am2 

ah4=ah2-an2 

ak= ( (ah4- (ah3*ag4/ag3) ) **2 . ) / (ag4* (ah3-ah4) ) 
as=ah3*ak 

ar=ak* (0 .25*ak- (ah3/ag3) ) 
ap=- . 25*ak- ( ,25*ak- (ah3/ag3) ) 
do 220 j=jdim+l,j2 
j j j=2* jdim-j 
xl (i, j) =xl (i, j j j) 

zl (i, j) =ap*xl (i, j) - ( (abs (ar* (xl (i, j) **2 . ) +as*xl (i, j) ) ) **0.5) 

220 continue 

do 230 j=jdim+l, j2 
x2 ( i , j ) =xl ( i , j ) +am2 
z2 (i, j ) =zl (i, j) +an2 
230 continue 

x2 (i, jdim) =x2 (i, jdim+1) 

z2 (i, jdim) = (z2 (i, jdim+1) +z2 (i, jdim-1) ) /2 . 
do 202 j=l, j2 

201 format (lx, 2 (i3, 2x) , lx, 4(fl4.4,lx) ) 

202 continue 
s (1)=0.0 

do 990 j=l, j2 
if (i .eq. 24) then 

if (j .It. j2) then 

dels= ( (x2 (i, j+1) -x2 (i, j) ) **2+(z2 (i, j+1) -z2 (i, j) ) **2) **.5 

s ( j+l)=s ( j) +dels 

endif 

if (j .eq. j2) then 

generates slopes and curvatures at any particular cross section 

hl=x2 (i, j-1) -x2 (i, j) 
h3=x2 (i, j-2) -x2 (i, j) 
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107 

990 

444 

200 


down2=hl*h3*(hl-h3) 

up2= (h3**2) *z2 (i, j-1) - (hl**2) *z2 (i, j-2) - ( (h3**2) - 
1 hl**2 ) *z2 (i, j ) 

slopel (i, j) =up2/down2 
denom2=hl*h3* (hl-h3) /2 

anumer2=h3*z2 (i, j-1) -hl*z2 (i, j-2) - (h3-hl) *z2 (i, j) 
curv (i, j) = (anumer2/denom2) / ( (1+slopel (i, j) **2) **1 .5) 
endif 

if (j .eq. 1) then 
hl=x2 (i, j) -x2 (i, j+1) 
h3=x2 (i, j) -x2 (i, j+2) 
downl=hl*h3* (h3-hl) 

upl— (h3**2) *z2 (i, j+1) - (hl**2 ) *z2 (i, j+2) - ( (h3**2) - 
1 hl**2) *z2 (i, j ) 

slopel (i/ j) =upl/downl 
denoml=hl*h3* (hl-h3) /2 

amimerl=h3*z2 (i, j+1) -hl*z2 (i, j+2) - (h3-hl) *z2 (i, j) 
curv(i, j ) = (anumerl/denoml) / ( (1+slopel (i, j) **2) **1.5) 
endif 

if ((j .ge. 2) .and. (j .It. j2) ) then 
h2=x2 (i, j) -x2 (i, j+1) 
hl=x2 (1, j-1) -x2 (i, j) 

slopel (i, j ) = (z2 (i, j+1) -z2 <i, j-1) ) / (x2 (i, j-1) -x2 (i, j+1) ) 
denom= ( (hi* (h2**2) ) + (h2* (hl**2) ) ) /2 

anumer= (hl*z2 (i, j+1) ) + (h2*z2 (i, j-1) )-( (hl+h2) *z2 (i, j) ) 
curv(i, j) = (anumer/denom) / ( (1+slopel (i, j) **2) **1 .5) 
write (13,*)hl,h2, -x2 (i, j) , z2 (i, j) , slopel (i, j ) , curvfi, j ) 
endif 

endif 

f ormat (2 (f 14 . 9, lx) ) 
continue 

if (i .eq. 24) then 

do 444 j=l, j2 

ratio=s ( j) /s ( j2) 

write (14, 107) ratio, slopel (i, j) 

write (15,107) ratio, curv (i, j) 

continue 

endif 

continue 

do 250 j=l,j2init 
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250 


z2 (1, j)=-8.589 
continue 
do 260 i=l, 32 
if (i . le. 22) then 

jdim= jdiminit 
j2= j2init 
else 

j dim= j dimint e r 
j2= j2inter 

endif 

do 255 j=l, j2 
if (j .It. ( jdim-1) ) then 

jj = j 

xx2 (i, j j) =x2 (i, j) 
zz2 (i, j j)=z2 (i, j ) 

endif 

if (j .eq. (jdim-1)) go to 255 
if (j .eq. jdim) then 
jj= jdim-1 
xx2 (i, j j)*=x2 <i, j) 
zz2 (if j j) =z2 (i, j) 

endif 

if (j .eq. (jdim+1)) go to 255 
if (j .gt. (jdim+1)) then 

j j-j-2 

xx2 (i f j j)=x2 (i f j ) 
zz2 (if j j ) =z2 (if j) 

endif 
255 continue 

260 continue 

do 263 i=23, 32 
do 264 j=l,j2init 
x2 (i f j ) =xx2 (i f j) 
z2 (i, j) =zz2 (if j) 

264 continue 

263 continue 

do 721 i=23, 44 
jextra=3 

topdels= (z2 (i f jdimf inal) -z2 (i f jdimfinal-1) ) / ( jextra+1) 
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topdelr= (x2 (i, jdimfinal) -x2 (i, jdimf inal-1) ) / ( jextra+1) 
botdelr= (x2 (i, jdimf inal+1) -x2 (i, jdimfinal) ) / (jextra+1) 
botdels= (z2 (i, jdimfinal+1) -z2 (i, jdimfinal) ) / (jextra+1) 
do 723 j=l, j2f inal 
if (j .le. ( jdimf inal-1) ) then 

j j=j 

xx2 (i f j j)=x2 (i, j) 
zz2 (if j j) =z2 (i, j) 

endif 

if (j .eq. jdimfinal) then 
j j= jdimf inal+jextra 
xx2 (i f j j)=x2 (i f j ) 
zz2 (i f j j)=z2 (i f j ) 

endif 

if (j .ge. (jdimfinal+1)) then 
j j-j+ (2* jextra) 
xx2 (i f j j)=x2 (if j) 
zz2 (if j j ) =z2 (i, j) 

endif 

723 continue 

do 726 jibba=l, jextra 
jj= jdimf inal-l+jibba 

xx2 (i f j j) «x2 (i f jdimf inal-1) + ( jibba*topdelr) 
zz2 (i f j j) =z2 (i, jdimfinal-1) + ( jibba*topdels) 

726 continue 

do 728 jibbal=l f jextra 
jj= jdimf inal+jextra+jibbal 

xx2 (i, j j) =x2 (i, jdimfinal) + ( jibbal*botdelr) 
zz2 (i f j j)=z2 (i, jdimfinal) +(jibbal*botdels) 

728 continue 

721 continue 

do 300 i=l f 44 
ylangley= (y (i) -42.5) *.17 
do 300 j=l, j2init-2 
yyy (if jf D = (y (i) -42.5) *.17 
300 continue 

do 500 i=l,44 

do 400 j=l f j2init-2 

xxx (i, j / l)=-xx2(if j) *.17 
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zzz (i, j, 1) =zz2 (i, j) * . 17 
400 continue 

500 continue 

do 7000 ii=l , iext 

yext (ii) =yahead*sin ( ( (ii-1) *pi) / (2* (iext-1) ) ) 

if (ii .eq. iext .or. ii .eq. (iext-1)) go to 7000 

do 7100 j=l, j2init-2 

assl=0 . 0 

ass2=0 . 0 

ass3=yext (ii) * . 17 

write (12, 530) as si, ass2, ass 3 

if ((j .eq. 1) .or. (j .eq. ( j2init-2) ) ) then 

write (16,530)assl,ass2,ass3 

endif 

7100 continue 
7000 continue 
c 

c writes input to grid generator ( GRIDTOOL ) 

c 

do 510 i=l, 22 

if ( (i .eq. 1) .or. (i .eq. 22) .or. (i .eq. 18)) then 
if (i .eq. 18) then 
dely=y (i) -y (1) 
endif 

if (i .eq. 22)then 
dely=y (i) -y (18) 
endif 

do 520 j=l,j2init-2 

assl=-xx2 (i, j) * . 17 

ass2= (zz2 (i, j) -zz2 (1, j) ) *.17 

ass3=(y(i) -42.5+100) *.17 

write (12,530) assl, ass2, ass 3 

if ((j .eq. 1) .or. (j .eq. ( j2init-2) ) ) 

1 write (16, 530) assl, ass2,ass3 
c if (i .eq. 18) then 

c write (8, *) assl, ass2 

c endif 

530 format (3f 15. 8) 

520 continue 

endif 
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510 continue 

do 550 i=23, 44 
dely=y (i) -y (i-1) 
c write (6, *) dely 

do 560 j=l, j2init-2 

assl=-xx2 (i, j) * . 17 

ass2= (zz2 (i f j) -zz2 (1, j) ) * .17 

ass3= (y (i) -42 . 5+100) * . 17 

write (12, 570) assl, ass2, ass3 

if ((j .eq. 1) .or. ( j . eq. ( j2init-2) ) ) 

write (16, 570) assl, ass2, ass 3 

endif 

570 format (3f 15. 8) 

560 continue 

550 continue 

go to 1000 
idiml=44 
jdiml= j2init-2 
kdiml=l 
c 

c write output in PLOT3D format 

c 

write(10,*) idiml , j diml , kdiml 

write (10, *) ( ( (xxx(i, j,k) ,i-l, idiml) , j=l, jdiml) , k= 
write (10, *) ( ( (yyy (i, j , k) , i=l, idiml) , j=l, jdiml) ,k= 
write (10, *) ( ( (zzz (i, j,k) ,i=l, idiml) , j=l, jdiml) ,k= 
1000 continue 

stop 
end 


l f kdiml) 
1, kdiml) 
1, kdiml) 
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APPENDIX C 


ERICKSON FOREBODY MATH MODEL 


The Erickson forebody math model was developed using the descriptions of the model contained in 
the Langley reports, a CAD drawing of the forebody and digitized data of the drawing generated by 
Nielsen Engineering Company. 


c Program efbmm.f (Erickson forebody math model) 


c Erickson/Hall chine forebody surface geometry definition 
c 

c W.H. Mason, December 1989 

c the main program is an example driver for the five subroutines 

c longitudinal lines routines: 

c xstn is input: ymhb, zupcl, and zlowcl are output 

c subroutine ymhbc (xstn, ymhb) - find the y max half breadth line 

c subroutine topcl (xstn, zupcl) - find the top center line 

c subroutine hotel (xstn, zlowcl) - find the bottom center line 

c 

c cross section lines routines: 


c 

c 

c 

c 

c 

c 


xstn, zupcl, zlowcl, ymhb and xs are input: zu and zl are output 

subroutine top (xstn, zupcl, ymhb, xs, zu) 
subroutine bot (zlowcl, ymhb, xs, zl) 

input /output is dimensional except xs, which is the % spanwise location 

xs must be positive, use symmetry to get lines for -xs 

the routines use spline subroutines from Conte and Deboor (included) 

dimension xpt (10,221) ,zpt(10,221) ,xstref (20) 


imax 

xstref (1) 
xstref (2) 
xstref (3) 
xstref (4) 
xstref (5) 
xstref (6) 
xstref (7) 
xstref (8) 


= 9 

= 0.0 

= 1.44 

= 4.312 

= 7.19 

= 10.12 
= 13.56 
- 17.25 
= 19.94 
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xstref (9) 


30.00 


do 10 istn = l,imax 

xstn = xstref (istn) 

call ymhbc (xstn, ymhb) 
call topcl (xstn, zupcl) 
call hotel (xstn, zlowcl) 

write ( 6, 80) istn, xstn, ymhb, zupcl, zlowcl 

do 10 i 1,101 

ii = 201 - (i-1) 

xs = 0.01*float (i-1) 

call top (xstn, zupcl, ymhb, xs, zu) 
call bot (zlowcl, ymhb, xs, zl) 
xpt(istn,i) = xs*ymhb 
zpt ( istn, i) = zu 
xpt(istn,ii) = xs*ymhb 
zpt (istn, ii) = zl 
10 continue 

do 20 i = 1,201 

write(6,100) (xpt ( j , i) , zpt ( j , i) , j = l,imax) 
20 continue 

80 format (2x, i4, f 10 .5, 3f 15 .5) 

100 format (18f 14 . 6) 
stop 
end 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


subroutine top (xstn, topcl, ymhb, xs, zu) 

W.H. Mason, December, 1989 

upper surface chine cross-section based on Bob Hall's routine 

xs is % semi-span, and zu is dimensional 

topcl - top centerline of body 

ymhb - max-half breadth of configuration 

xs - spanwise location where z value is desired 

zu - output value of body elevation at input xs 

dimension xi (11) , c (4, 11) 

define cross section for xstn up to 7.19 


Pi 

= 3.1415926 

zu 

= 0.0 
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if (xstn . le. 0.0) return 


c 

c 

c 

c 


theta = 12 . 0*pi/180 . 0 


cc 

tantht 

a4 

a3 

a2 


= topcl/ymhb 
= tan (theta) 

= 2.*tantht - 3.*cc 
= 5.*tantht - 8.*cc 
= 3.*tantht - 6.*cc 


zul 


(a4*xs**4 - a3*xs**3 + a2*xs**2 + cc) * 


define cross section for xstn = 19.94 and beyond 


if(xs .le. 0.40) then 


zu2 = sqrt (1.0-2. 6*xs**2) *topcl 
go to 10 

end if 


nin 

= 

4 

xi(l) 

= 

.4 

xi(2) 


0.62594 

xi (3) 


0.90121 

xi (4) 

= 

1.0 

c(l,l) 

= 

0.764199 

c (1, 2) 

= 

0.34556 

c (1, 3) 

= 

0.022085 

c(l,4) 

= 

0.0 

c(2,l) 

= 

-1.3609 

c (2, nin) 

= 

-0.175 


n = nin - 1 

call spline ( n , xi, c) 
call calccf (n,xi, c) 

zu2 = pcubic (xs,n,xi, c) *topcl 

10 continue 

if (xstn .le. 7.19) zu = zul 
if (xstn .ge. 19.94) zu = zu2 

if (xstn .gt. 7.19 .and. xstn .It. 19.94) then 

xn = (xstn - 7.19)/ (19 . 94 - 7.19) 
cn = f (xn) 

zu = (1. - cn) *zul + cn*zu2 


ymhb 
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end if 


return 

end 

function f (xn) 
c blending function 

y = 5 . 0*xn**2 - 6.0*xn**3 

if (xn .gt. 0.5) y = 1 . - (5 . 0* (1 . -xn) **2-6 . * (1 . -xn) **3) 

f = y 

return 

end 

subroutine bot (botcl, ymhb, xs, zl) 
c 

c W. H. Mason, December, 1989 

c 

c lower surface chine cross-section based on Bob Hall's routine 

c 

c xs is % semi-span, and zl is dimensional 
c botcl - bottom centerline of body 

c ymhb - max-half breadth of configuration at this station 

c 

c xs - spanwise location where y value is desired 

c zl - output value of body elevation at input xs 

dimension xi (60) , c (4, 60) 

pi = 3.1415926 

zl = 0.0 

if (ymhb .eq. 0.0) return 

thetal = -2 . 0*pi/180 . 0 

b = botcl/ymhb 

tantht « tan (thetal) 

a 4 = 3 . 0*b - 2 . *tantht 

a3 = 8.0*b - 5.*tantht 

a2 = 6 . 0 *b - 3 . *t antht 

c 

zlbase = (a4*xs**4 - a3*xs**3 + a2*xs**2 - b) * ymhb 

zmod = 0.0 

nin - 13 

xi ( 1 ) = 0.000 


62 



xi<2) 

= 

0.020 

xi (3) 

= 

0.120 

xi(4) 

= 

0.230 

xi (5) 

= 

0.340 

xi (6) 

= 

0.450 

xi (7) 

= 

0.540 

xi (8) 

= 

0.590 

xi (9) 

= 

0.630 

xi (10) 

= 

0.740 

xi(ll) 

= 

0.850 

xi (12) 

= 

0.940 

xi (13) 

= 

1.000 

c(l,l) 


0.000000 

C dr 2) 

= 

0.002180 

c (1,3) 


0.050257 

c (1, 4) 

= 

0.152545 

c (1, 5) 


0.249569 

c (1, 6) 

= 

0.300617 

c (1, 7) 

= 

0.270427 

c (1, 8) 

= 

0.202228 

c (1, 9) 

= 

0.135704 

c (1, 10) 

= 

0.050513 

c(l,ll) 


0.017489 

c (1, 12) 

= 

0.003740 

c (1, 13) 


0.000000 

c (2, 1) 

= 

o 

o 

c (2, nin) 

= 

0.0 

n 

— 

nin-1 


call spline (n, xi, c) 
call calccf (n, xi f c) 


bot719 = -1.6359 

c the 1.76794 factor is the ymhb at x = 7.19 

zmod = pcubic (xs, n, xi, c) * (botcl/bot719) *1.76794 

zl = -(zlbase + zmod) 

return 

end 


subroutine ymhb c (xstn, ymhb) 
c 

c y max half breadth definition for the Erickson chine forebody 
c 

c xstn is input , ymhb is the output 
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c W.H. Mason, final version May 16, 1990 
dimension xi (20) , c (4, 20) 
if (xstn . le. 0.0) ymhb = 0.0 

c tangent ogive portion 

if (xstn .gt. 0.0 .and. xstn .le. 8.0) then 


rO 

xlref 

rOxl 

deltana 

deltan 

xnrO 

xnxl 

a 

b 

xn 

x 

xl 

rxl 

ymhb 

end if 


- 4.35 
= 27.35 
= rO/xlref 
= 27.5 

= deltana*3. 1415926/180.0 
= sin (deltan) / (1 . 0 - cos (deltan) ) 

= xnr0*r0xl 

= r0xl*cos (deltan) / (1 . 0 - cos (deltan) ) 
= 2 . *r0xl*sin (deltan) / (1 . 0 - cos (deltan) ) 

= xnxl*xlref 

= xstn 
= x/xlref 

= sqrt(a**2 + b*xl - xl**2) - a 
= rxl * xlref 


c first spline pfortion 

if (xstn .gt. 8.0 .and. xstn .It. 19.0) then 


nin 


3 

xi(l) 

= 

8.0 

xi (2 ) 

ss 

11.00 

xi (3) 

= 

19.00 

c(l,l) 

= 

3.08787 

c (1/ 2) 

= 

3.717 

c (1, 3) 

= 

4.35 

c (2, 1) 

= 

0.26256 

c (2, nin) 

= 

0.0 


n = nin - 1 

call spline (n, xi, c) 
call calccf (n,xi,c) 


64 



ymhb 


= pcubic (xstn, n, xi, c) 


end if 

c straight line portion 

if (xstn .gt. 19.0) ymhb = 4.35 

return 

end 


subroutine topcl (xstn, zupcl) 
c 

c top centerline definition for the Erickson chine forebody 

c 

c xstn is input, zupcl is output 

c W.H. Mason, final version May 165, 1990 

dimension xi (20) , c (4, 20) 

if (xstn .le. 0.0) zupcl = 0.0 

c tangent ogive portion 

if (xstn .gt. 0.0 .and. xstn .le. 7.5) then 

r0 = 2.616 

xlref = 23.000 

rOxl = rO /xlref 

deltana = 19.500 

deltan = deltana*3 . 1415926/180 . 0 

xnrO = sin (deltan) / (1 . 0 - cos (deltan) ) 

xnxl = xnr0*r0xl 

a = r0xl*cos (deltan) / (1 . 0 - cos (deltan) ) 

b = 2 . *r0xl*sin (deltan) / (1 . 0 - cos (deltan) ) 

xn = xnxl*xlref 

x = xstn 

xl = x/xlref 

rxl = sqrt (a**2 + b*xl - xl**2) - a 

zupcl = rxl * xlref 

end if 

c first spline ptortion 
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if (xstn .gt. 7.5 .and. xstn .It. 14.94) then 


nin 

= 

3 

xi(l) 

=5 

7.5 

xi(2) 

= 

10.12 

xi (3) 

= 

14.94 

c(l,l> 

= 

1.9571 

c (1, 2) 

= 

2.35 

c (1, 3) 

= 

2.58 

c(2,l) 

= 

.17184 

c (2, nin) 

= 

0.0064 


n = nin - 1 

call spline (n, xi, c) 
call calccf (n, xi, c) 

zupcl = pcubic (xstn, n, xi, c) 

end if 

c first straight line portion 

if (xstn . ge . 14.94 .and. xstn .le. 17.25) then 

xref = 14.94 

zref = 2.58 

dzref = 0.0064935 

zupcl = zref + (xstn - xref) *dzref 

end if 

c second spline section 

if (xstn .gt. 17.25 .and. xstn .It. 19.94) then 


nin 

= 2 

xi (1) 

= 17.25 

xi (2) 

= 19.94 

c(l,l) 

= 2.595 

c (1/ 2) 

- 2.630 

0(2,1) 

= 0.0064935 

c (2, nin) 

= 0.0294118 


n = nin - 1 
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call spline (n, xi,c) 
call calccf (n, xi, c) 

zupcl = pcubic(xstn, n, xi,c) 

end if 


c 


second straight section 


if (xstn .ge. 19.94 .and. xstn .le. 23.0) then 


xref 

zref 

dzref 


19.94 

2.63 

0.0294118 


zupcl = zref + (xstn - xref) *dzref 
end if 


c third and final spline section 

if (xstn .gt. 23.0 .and. xstn .le. 30.0) then 


nin 

= 

2 

xi (1) 


23.0 

xi (2) 

= 

30.0 

c(l,l) 

= 

2.72 

c (1/2) 

= 

2.75 

c (2, 1) 


0.0294118 

c (2, nin) 

= 

0.0 

n 

= 

nin - 1 


call spline (n, xi, c) 
call calccf (n, xi, c) 

zupcl = pcubic (xstn, n, xi, c) 

end if 

if (xstn .gt. 30.0) zupcl = 2.75 

return 

end 

subroutine botcl (xstn, zlowcl) 

bottom centerline definition for the Erickson chine forebody 
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c xstn is input, zlowcl is output 

c 

c W.H. Mason, final version May 16, 1990 

dimension xi (20) , c (4, 20) 

if (xstn .le. 0.0) zlowcl = 0.0 

c tangent ogive portion 

if (xstn .gt. 0.0 .and. xstn .le. 6.0) then 

r0 = 3.50 

xlref = 27.35 

rOxl = rO/xlref 

deltana = 15.0 

deltan - deltana*3 . 1415926/180 . 0 

xnrO = sin (deltan) / (1 . 0 - cos (deltan) ) 

xnxl « xnr0*r0xl 

a = r0xl*cos (deltan) /( 1 . 0 - cos (deltan) ) 

b = 2 . *r0xl*sin (deltan) / (1 . 0 - cos (deltan)) 

xn = xnxl*xlref 

x = xstn 

xl = x/xlref 

rxl = sqrt(a**2 + b*xl - xl**2) - a 

zlowcl — -rxl * xlref 

end if 

c first spline p[ortion 

if (xstn .gt. 6.0 .and. xstn .It. 19.94) then 


nin 

= 

3 

xi (1) 


6.0 

xi (2) 

= 

10.12 

xi (3) 

= 

19.94 

c(l,l) 

= 

-1.41617 

eg 

H 

0 


-2.025 

c (1, 3) 

= 

-2.59 

c (2, 1) 

= 

-0.20456 

c (2, nin) 

= 

-0.052288 


n = nin - 1 

call spline (n,xi, c) 
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call calccf (n, xi, c) 


c 


zlowcl = pcubic (xstn, n, xi,c) 
end if 


first straight line portion 
if (xstn .ge. 19.94 .and. xstn .le. 23.0) then 


xref = 19.94 

zref = -2.59 

dzref = -0.052288 


zlowcl = zref + (xstn - xref) *dzref 
end if 


c second spline section 

if (xstn .gt. 23.0 .and. xstn .It. 30.0) then 


nin 

= 

2 

xi (1) 

= 

23.0 

xi (2) 

= 

30.0 

c(l,l) 

= 

-2.75 

c(l,2) 

= 

-2.85 

c (2, 1) 

= 

-0.052288 

c (2, nin) 

= 

0.0 

n 

= 

nin - 1 


call spline (n, xi, c) 
call calccf (n,xi,c) 

zlowcl = pcubic (xstn, n,xi, c) 

end if 

c second straight section 

if (xstn .ge. 30.0) then 

zlowcl = -2.85 
end if 

return 

end 
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SUBROUTINE SPLINE (N, XI, C) 

C 

C FROM CONTE AND DEBOOR 

C 

DIMENSION XI (201) , C (4,201) ,D(201) ,DIAG(3201) 
DATA DIAG(l) ,D(1) /l . 0D0, 0 . 0D0/ 

NP1=N+1 
DO 10 M=2,NP1 
D(M)=XI (M) -XI (M-l) 

10 DIAG (M) = (C (1,M) -C (1,M-1) ) /D (M) 

DO 20 M=2,N 

C(2,M)=3.0D0*(D(M) *DIAG (M+l) +D (M+l ) *DIAG (M) ) 

20 DIAG (M) =2 . 0D0* (D (M) +D (M+l) ) 

DO 30 M=2 , N 

G=-D (M+l) /DIAG (M-l) 

DIAG (M) =DIAG (M) +G*D (M-l) 

30 C (2,M) =C (2,M) +G*C (2, M-l) 

NJ=NP1 
DO 40 M=2, N 
NJ=NJ-1 

40 C (2,NJ) = (C (2,NJ) -D(NJ) *C (2,NJ+1) ) /DIAG (NJ) 
RETURN 
END 

SUBROUTINE CALCCF (N, XI, C) 

DIMENSION XI (201) ,C(4,201) 

DO 10 1=1, N 
DX=XI (1+1) —XI (I) 

DIVDF1=(C(1,I+1)-C(1,I))/DX 
DIVDF3=C (2, I) +C (2,1+1) -2 . 0D0*DIVDF1 
C (3, I) = (DIVDF1-C (2, I) -DIVDF3) /DX 
10 C (4, I) =DIVDF3/DX/DX 
RETURN 
END 

FUNCTION PCUBIC (XBAR, N, XI, C) 

DIMENSION XI (201) ,C(4,201) 

1 = 1 

DX=XB AR-X I (I) 

IF (DX) 10,30,20 
10 IF (I.EQ.l) GO TO 30 
1 = 1-1 

DX=XBAR-XI (I) 

IF (DX) 10,30,30 

19 1=1+1 
DX=DDX 

20 IF(I.EQ.N) GO TO 30 
DDX=XBAR-XI (1+1) 

IF(DDX) 30,19,19 

30 PCUBIC=C (1, 1) +DX* (C (2, 1)+DX* (C(3, 1)+DX*C (4, 1) ) ) 
RETURN 
END 
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TANGENT OGIVE 3 1 x 69 x 24 grid 



a = 20° ,P = 0° 

a = 0° , (3 = 20° 

Cx 

-0.02004 

-0.02004 

c y 

0.2737E-3 

0.1842 

Cz 

0.1845 

-0.3878E-4 


-1.595 

0.2988E-3 

c n 

-0.4143E-2 

1.593 

Iterations 

1950 

2260 


Table 1. Comparison of Forces and Moments showing equivalence between angle-of-attack 
with no sideslip and sideslip with no angle-of-attack for an axisymmetric body 



a - 50 ° - 60 ° 



Figure 1. Typical C nR characteristics of an advanced fighter 




Figure 3. Schematic of a typical forebody flowfield characteristics 



Figure 4. Goal of forebody design effort 
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(b) PLAN AND PROFILE VIEWS 


Figure 5. Samples of shapes employed in fighter forebody design 
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Figure 6. Tangent ogive crude grid Euler solution comparison with experimental data on 
windward plane. 
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Figure 7(a). Tangent ogive crude grid Euler solution comparison with 
experimental data(x/d = 0.5) 
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Figure 7(b). Tangent ogive crude grid Euler solution comparison with 
experimental data (x/d = 2.0) 
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Figure 7(c). Tangent ogive crude grid Euler solution comparison with experimental data 
(x/d = 3.5) 



theto 

Figure 7(d). Tangent ogive crude grid Euler solution comparison vrith experimental dam 
(x/d = 6.0) 


79 




.5 



theta 


Figure 8(a). Tangent ogive : Demonstration of a / f} results at x/d = 0.062 (Re=2*10^) 




Figure 8(b). Tangent ogive : Demonstration of a / /? results at x/d = 1 1.02 (Re=2* 10^) 



Figure 8(c). Tangent ogive : Demonstration of a/p results at x/d = 24.0 (Re=2*10 5 ) 
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Figure 9(a). F-5A geometry ( from Tom Heglund of Northrop) 
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Figure 9(b). F5-A geometry: FS 49.70 in.(forward station) 
(from Tom Heglund of Northrop) 



Figure 9(c). F5-A geometry: FS 130.35 in.(just prior to.sft^ight'Jiiur side beginning) 
(from Tom Heglund of Northrop) 
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Figure 9(e). F5-A geometry: FS 194.00 in.(end station) 
(from Tom Heglund of Northrop) 
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Figure 10(a). F5-A wind tunnel model being measured to determine the 
shape. 
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Figure 10(b). F5-A wind tunnel model being measured to determine the 
shape (close-up). 


86 


-0.5 
- 1.0 
z > in -1.5 
- 2.0 
-2.5 

-1.5 -1 -0.5 0 0.5 1 1.5 

y, in 

Figure 1 1(a). Comparison of F-5A geometry: Station 3.375 in from nose 
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Figure 1 1 (b). Comparison of F-5A geometry: Station 6.25 in from nose 
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• W.T. model inspection contour 

o F-5A Theoretical surface 



Figure 11(c). Comparison of F-5A geometry: Station 9.1875 in from nose 


a F-5A theoretical surface 

• W.T. model inspection contour 



y, in 

Figure 11(d). Comparison of F-5A geometry: Station 15.3125 in from nose 
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• W.T. model inspection contour 

□ F-5 A Theoretical surface 



y, in 

Figure 11(e). Comparison of F-5 A geometry: Station 18.0625 in from nose 


• W.T. model inspection contour 

n F-5 A theoretical surface 



y, in 

Figure 1 1(f). Comparison of F-5A geometry: Station 22.5 in from nose 
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Figure 12. Wireframe model of F-5A surface definition 
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(a) cross sectional grid 



(b) closeup near surface 

Figure 13. F5-A cross sectional grid and closeup near surface at x = 9.061 in. 
( model scale ) 
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(b) closeup near surface 

Figure 14. F5-A cross sectional grid and closeup near surface at x = 15.44 in. 
( model scale ) 
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(a) cross sectional grid 
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Figure 16. F-5A longitudinal pattern of grid (PL0T3D) 
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Figure 20(a). F5-A inviscid surface pressure distribution at x = 14.02 
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(b)x = 29.61 


licure 22. F5-A inviscid vs turbulent surface pressure distribution at 
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Figure 24. F5-A inviscid pressure contours at x = 14.025 in.( a = 40° and (5 = 5°) 



Figure 25. F5-A turbulent pressure contours at x = 14.025 in.( a = 40° and (3 = 5°) 
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Figure 26. F5-A inviscid stagnation pressure at x = 14.025 in.( a = 40° and = 5°) 



Figure 27. F5-A turbulent stagnation pressure at x = 14.025 in.( a = 40° and /? = 5°) 
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Figure 28. Details of the Erickson forebody wind tunnel model ( ref. 7 ) 
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Figure 31. Erickson forebody computational model (PLOT3D) 
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(b) closeup near surface 

Figure 32. Erickson forebody cross sectional grid at x = 30.00 in 
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Figure 35(c). Erickson forebody inviscid surface pressure distribution at x 

( M = 0.2) 


= 19.94 in. 
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Figure 37(a). Erickson forebody inviscid vs turbulent surface pressures at x = 7.19 in. 
( M = 0.2, Re c = 1.02*10 6 ) 
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Figure 37(b). Erickson forebody inviscid vs turbulent surface pressures 
( M = 0.2, Re c = 1.02*10 6 ) 


atx = 13.56 in 
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Figure 37(c). Erickson forebody inviscid vs turbulent surface pressures at x = 19.94 in. 
( M = 0.2, Re c = 1.02*10 6 ) 
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Figure 38(a). Erickson forebody invisc 


( M = 0.2, Re c = 1.02*10 








Figure 38(c). Erickson forebody inviscid 
( M = 0.2, Re c = 1.02*10 6 ) 
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Figure 39. Erickson forebody with isolated wing VLM model 
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Figure 40. Wing induced flowfield at Erickson forebody pressure stations. 
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Figure 41(a). Erickson forebody inviscid surface pressures at x = 7.19 in. 
(a =32.2° and ft =0°) 
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Figure 41(b). Erickson forebody inviscid surface pressures at x = 13.56 in. 
(a = 32.2° and B = 0°) 
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Figure 41 (c). Erickson forebody inviscid surface pressures at x = 19.94 in. 
(a =32.2° and 0=0°) 
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Figure 42. Erickson forebody inviscid surface oil flow pattern (a = 30° and /? = 0°) 



Figure 43. Erickson forebody turbulent surface oil flow pattern (a = 30° and p = 0°) 
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Figure 44. Erickson forebody inviscid pressure contours at x = 13.56 in. 
(a =30° and /J = 0°) 



Figure 45. Erickson forebody turbulent pressure contours at x = 13.56 in. 
(a =30° and /? = 0°) 
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Figure 46. Erickson forebody inviscid stagnation pressure at x = 13.56 in. 
(a =30° and 0=0°) 



Figure 47. Erickson forebody turbulent stagnation pressure at x=13.56in. 
(a =30° and/3 =0°) 
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Figure 48. Erickson forebody inviscid surface flow pattern ( a = 30° and (3 = 10° ) 
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Figure 49. Erickson fore body turbulent surface flow pattern ( a = 30° and (3 = 10° ) 
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Figure 50. Erickson forebody inviscid pressure contours at x = 13.56 in. 
(a =30° and /? = 10°) 


Figure 51. Erickson forebody turbulent pressure contours at x = 13.56 in. 
(a = 30° and ft = 10°) 
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Figure 52. Erickson forebody inviscid stagnation pressure at x = 13.56 in. 
(a =30° and/3 = 10°) 



Figure 53. Erickson forebody turbulent stagnation pressure at x=13.56in. 
(a = 30° and p = 10°) 
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Figure 55. Parametric forebody shapes when the straight sidewall is selected. 
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Figure 56. Examples of the parametric forebody for the entire cross section 
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